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Introduction 

Since the launching of the first Sputnik in 1957, the use of near- 
Earth satellites for scientific, military, and commercial applications has 
progressed far beyond the expectations of even the most visionary of the 
early satellite pioneers. As a result, many sophisticated numerical and 
mathematical techniques have evolved to permit the precise prediction 
of the position of a satellite in time and space, and the Earth s external 
gravitational field has been determined to high accuracy. 

However, for most mission planning purposes, and even in fact for 
many operational applications, ultra-high precision is neither warranted 
nor desired, as the computational expense generally increases quite 
rapidly with increased computational precision. Consequently, there 
is still a firm need for simple computational algorithms which produce 
results of reasonable accuracy at moderate expense. 

Anyone who works in a specialized scientific discipline for any length 
of time naturally acquires over the years a computational methodol- 
ogy, consisting of large numbers of computational algorithms of many 
degrees of complexity and accuracy, other computational procedures, 
approximate methods, etc., with which he feels most comfortable and 
which have proven to be useful and accurate, even though they may 
not always provide the most direct way of making a particular calcu- 
lation. The present text is such a compilation of procedures in orbital 
mechanics and spherical astronomy which the author has collected over 
about a 30-year period. Many of the formulas presented herein were 
included in an attempt to make the text as self-contained as possible, 
without requiring the reader to consult almanacs, ephemerides, or other 
reference material. Nothing but the material presented in this report is 
necessary to carry out any of the computations described, thus, these 
procedures are ideal for programming on computers. 

It is not the author’s intention to derive or otherwise develop all 
the equations and/or formulas presented herein — this task is most ade- 
quately covered ip practically any good textbook on Celestial Mechan- 
ics or Spherical Astronomy. (See, for example, any of the general texts 
cited.) Any equations, however, which are not widely circulated in the 
literature, are either specifically referenced or developed to the point 
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where a heuristic argument will suffice to convince the reader that the 
equation does what it is claimed to do and what its limitations are. All 
the comments, opinions, or other commentary are, of course, the au- 
thor’s own, and perhaps serve to do nothing but illustrate the author’s 
own prejudices. 

The problem of presenting a completely unambiguous and unique 
set of symbols is one that has plagued writers of scientific (and pseudo- 
scientific) documents for many years, and the present author offers no 
relief in this regard. In several instances, the same symbol may stand 
for two or more completely unrelated quantities (e.g., a describes a 
coordinate axis, also stands for the component of a vector along this 
axis, and is used as the symbol for the zenith distance angle). The 
philosophy adopted here is that, when discussing a particular area or 
discipline, the symbols most commonly used by most writers in that 
area are used. These are defined in the text as encountered and the 
use of multiply defined symbols should cause no problem as they are 
generally used only in the few pages in the neighborhood of where they 
are defined. 

The layout of the text is as follows: the coordinate systems in 
which the remainder of the text is developed and the systems of time 
measurement common in the computation of satellite ephemerides are 
presented in chapter 1. The mathematical description of the external 
gravity field of the Earth is presented and the restriction in the text to 
the use of zonal harmonics is defended in chapter 2. Chapter 3 defines 
the author s preference for the set of orbital elements used to describe 
and advance in time the position and velocity of the spacecraft, while 
chapter 4 introduces the Cartesian form of the equations of motion. The 
Lagrange Planetary Equations are not used, but their use in theoretical 
developments is mentioned briefly, and some of these theoretical results 
are used herein. Numerical comparisons between sundry numerical 
procedures are presented. Chapter 5 presents numerical algorithms for 
accurately determining the right ascension and declination of the Sun. 
The word accurately here refers to the precision of computing a given 
parameter relative to that usually required by most near- Earth satellite 
applications and not generally to the precision available for, say, some 
astronomical applications. 

The text concludes with a few very brief remarks on the application 
of the present methods to other Earth satellite problems and to the 
computation of the orbits of the planets and the Moon. 



Chapter 1 

Coordinate Systems and Time 

The primary coordinate system used in the present text is a quasi- 
inertial system, defined by the Earth’s equator and the apparent orbit 
of the Sun around the Earth with the origin at the center of mass of the 
Earth. (See fig. 1-1.) The intersection of the Sun’s orbital plane, the 
ecliptic , with the Earth’s equatorial plane defines a line, called the line 
of nodes. The direction defined by the center of the Earth and the node 
at which the Sun appears to cross the equator from south to north is 
called the ascending node , the first point of Aries , or more usually, the 
vernal equinox qf and defines the direction of the x-axis. The rotational 
axis of the Earth defines the 2 -axis, and the t/-axis is located in the 
equatorial plane in such a way that the xyz coordinate system is a 
right-handed one. The ecliptic plane is inclined to the equator at about 
23.44°, the obliquity of the ecliptic , and can be accurately computed 
from equations presented later in this chapter. 


z 



Figure 1-1. Definition of quasi-inertial x-axis by intersection of equatorial and 
ecliptic planes . 
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If the Earth and Sun were both perfectly homogeneous spheres so 
that their mutual gravitational attraction followed an inverse-square 
law and acted along a line joining their centers, if there were neither 
Moon nor planets, and if the Sun were either stationary in the universe 
or moved through it on a straight line at constant speed, then the 
coordinate system described in the previous paragraph would indeed be 
an inertial system. However, none of these conditions hold. The Earth’s 
figure, both geometrically and dynamically, is closely approximated not 
by a sphere but by an oblate spheroid — an ellipse rotated about its 
minor or shorter axis the equatorial radius being about 21 km larger 
than the polar radius (6378.160 km and 6356.775 km, respectively). As 
shown in figure 1-1, the ecliptic plane is inclined to the equator at about 
23.44°, and thus, the direction of the gravity resultant of the Earth-Sun 
does not pass through the center of mass of the Earth but instead passes 
through a point on the Sun side of the line joining the centers of the 
Earth and Sun. This produces a gravitational couple, or torque, on the 
Earth, and since the Earth is spinning, this torque causes the Earth to 
wobble or precess (fig. l-2(a)) and nutate (fig. l-2(b)). Likewise, the 
Earth has a rather large Moon which also orbits in a plane inclined 
to the equator and which also produces a gravitational torque. This 
combined lunisolar perturbation causes the Earth to precess and nutate, 
much as a spinning top does, and causes the 2-axis to rotate about the 
ecliptic pole K in a clockwise direction as seen from the North Pole. 
This, in turn, results in the vernal equinox moving clockwise in the 
plane of the ecliptic at the rate of about 50 arc-sec/yr, the precession of 
the equinoxes. This measured value also includes a small term, called 
the planetary precession, due to the other planets, which causes the 
plane of the ecliptic, assumed to be fixed when computing the lunisolar 
precession, to wobble slightly. The combined lunisolar and planetary 
precessions are called the general precession. 


(a) Precession only. ( b ) Precession plus nutation. 

Figure 1-2. Sketch illustrating precession only and precession plus nutation on 
motion of Earth’s pole. Motion of 1 is called precession of equinox; P is north 
pole of Earth; K is north pole of ecliptic; P moves around K. 
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Chapter 1 


When referring to astronomical coordinates, one generally fixes the 
vernal equinox (x-axis) and the equator in one of two ways — either by 
referring to their positions at a specific date and then requiring that 
these positions be frozen for the duration of the computational period, 
or by appending the words “of date” to the positions of the vernal 
equinox and equator, which means that these directions are slowly 
changing while the computations are being made and refer to their 
locations at the current time. If only the precessional effect is being 
taken into account (usual for much Earth satellite work) one refers to 
mean-of-date coordinates. If nutation is also included, then one refers to 
true-of-date coordinates. In the present work, mean-of-date coordinates 
are used, except where specifically noted, and are assumed synonymous 
with “inertial” coordinates. 1 

Celestial objects are located in this coordinate system by the use of 
two angles — right ascension and declination analogous to the more 
familiar longitude and latitude, respectively (fig. 1-3). The right 
ascension a is measured from the vernal equinox along the equator, 
positive east, or counterclockwise as seen from the North Pole of 
the Earth. The right ascension is frequently measured in time units 
(practically always in astronomy), but for computational purposes, 
degrees are used herein. The declination S is measured along a 
“meridian,” north and south from the equator positive north, and in 
the present work, since we will consider the Earth to be geometrically 
a sphere but gravitationally an oblate spheroid, the declination is 
numerically identical with the latitude. 

If R is the distance from the center of the Earth to the celestial 
object, then a vector to the body can be written in the right ascension- 
declination system as 


r = R 


cos 

cos 


<5 cos a 

6 sin a 

sin <5 


(M*) 


1 In chapter 5, which discusses the determination of the position of the Sun, 
proper allowance for the motion of the equinox is made in the nonlinear terms of 
the sundry time series used in the computations. Thus, the Sun’s position in the 
mean-of-date coordinates is accurately determined. No such allowance is made in 
the earlier chapters which discuss the position computation of an Earth satellite 
these coordinate systems are assumed to be truly inertial. The only orbital element 
that is affected by this assumption is the position angle of the line of modes (defined 
in chapter 4). A small linear term could be subtracted from the motion of this 
parameter computed in chapter 4, but the error induced in reflecting this effect is 
negligible compared with the inherent accuracy of these simplified equations. 
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Figure 1-3. Definitions of right ascension a and declination 6 of celestial body. 


from which the Cartesian coordinates of the body, in our inertial system, 
can be computed. Also, as will frequently be done, if the Cartesian 
coordinates are known, equation (1-la) can be used to compute the 
right ascension and declination. (Frequently, the distance of the body 
from the Earth center is either unknown or one has only a unit vector 
defining the direction of the body. Equation (1-la) can still be used to 
etermine the right ascension and declination simply by assuming that 
|r| and it = 1.) 

Measurements made on the Earth’s surface are most conveniently 
referred to a set of coordinate axes fixed with respect to the rotating 
Earth. The ar-axis is defined by the intersection of the meridian pass- 
ing through the Greenwich Observatory in England— known, cleverly 
enough, as the Greenwich meridian— with the equator. The rotating z- 
axis is again defined by the Earth’s spin axis, and the y-axis completes 
the right-handed system. Objects are located on the Earth by the fa- 
miliar latitude-longitude spherical coordinates. Latitude is measured 
positive northward from the equator, along a meridian. Longitude is 
measured in the equatorial plane, positive eastward from the Green- 
wich meridian. The point labeled P in figure 1-4 has latitude %p c and 
longitude A. The Greenwich meridian is also labeled, as is the position 
of the nonrotating rr-axis, which is labeled *, r Assuming a spherical 
Earth of radius R E} a vector to the surface point P can be written in 
the rotating coordinates 


r 


6 


r E = Re 


COS 

cos %p c 
sin 


cos 

sin 


A 

A 


(Mb) 
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Figure 1-4. Location of Greenwich meridian and definition of geocentric latitude 
t l> c and longitude Xg- Greenwich sidereal time is 0 g . 

In order to relate the coordinates of P in the Earth-fixed axes 
(eq. (2-lb)) to the inertial system (eq. (2-la)), the location of the 
Greenwich meridian at any time relative to the vernal equinox Og is 
needed. Then, in the inertial system J~ 



' x' 


’ COS Og 

— sin 0 g 

O' 


Re 

C08 

t \>c 

cos 

(Og + A) 

r = 

y 

= 

sin 0 g 

COS Og 

0 

r E — 

Re 

cos 


sin 



z 


0 

0 

1. 



Re 

sin 

li>C 

. 


The latitude t /> c , as determined by equation (1-2), is called the 
“geocentric” latitude of P (fig. l-5(a)) and is measured from the 
equatorial plane along a meridian to the line joining P and the origin 
of coordinates. 

The “geodetic” latitude xp g of figure 1-5 (b) is another measurement 
of latitude, usually but not always confined to measurements made on 
the Earth’s surface and is the angle between the equatorial plane and 
the normal to the geoidal surface which passes through P, line PAC. 
The altitude of a spacecraft at P is measured along the direction of 
PA, the local vertical or plumb bob direction, rather than along PB. In 
orbit mechanics, the point A is called the subsatellite point. 

If the point P were on the surface, as in figure 1-5 (a), there is a 
simple relation between the geocentric and geodetic latitudes, namely, 

2 

tan rp g = tan rp c (1*3) 
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(a) Point P on surface. 



(b) Spatial point P above surface. 

Figure 1-5 Sketch of geometry and definition of geodetic latitude +. and 
geocentric latitude ip c of point on and above surface of oblate spheroid. 
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Chapter 1 

where a and b are the equatorial and polar radii of the Earth, 6378.160 
and 6356.775 km, respectively. 

In the more general case shown in figure l-5(b), no such simple 
relation exists, neither for the latitudes nor for the height h min , and 
more complex methods must be resorted to. Escobal (1965) presents 
an iteration technique. The author’s own solution is presented here. 

In figure l-5(b), let Q be any point on the surface of the ellipse cross 
section. The distance between P and Q is then 


h 2 — (x - Xp ) 2 + (z — z p ) 2 


(1-4) 


The problem is to minimize h subject to the constraint that the point 
Q (x, z) lies on the ellipse 


x 2 * 2 _ , 

“o + To - 1 


b 2 


(1-5) 


This is a straightforward minimization problem requiring the use of the 
Lagrange multiplier technique. We want to form the function 



where g is the Lagrange multiplier. The three equations 


d$ , v 2 gx 

— = 2 (x - x p ) + = 0 

ox or 


f 


<94> , x 2 gz 

— = 2 (z - zp) + -p- = 0 


and equation (1-5) are sufficient to determine the three unknowns, 
x, z, and g. After a fair amount of elementary algebra we arrive at the 
relation 

F(0 = A 4 £ 4 + A 3 £ 3 + A 2 ? + Ai f + Aq (1-6) 
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where 

a 4 = q(i -r ) 2 
a 3 = 2<*r(i - r) 

A l 2 = ar 2 -(l-r ) 2 + /3 
Ai = -2r(i -r) 

A 0 = r 2 

*=(?) 2 

r =(?) 2 

u ~ X J Xp , anc * ^ = z l z v Equation (1-6) is a quartic equation 
which can be solved exactly using the classical Descartes’ method (See 
for example, Escobal 1965.) However, since the Earth is very nearly 
spherical, the author has found that the Newton-Raphson iteration 
technique generally converges in one or two iterations and prefers that 
method. In either case, once £ is found, then 

1 

V = t r 

from which 

x = $x p 1 

z = r l z p j 

give the coordinates for point A, the subsatellite point. 

"min 1S then found by substituting these coordinate 
equation (1-4). The geocentric latitude is found from 

ip c = tan -1 1 (i- 9 ) 

and the geodetic latitude from equation (1-3). 

Sidereal Time 

The angle 6 g in equation (1-2) is a very important quantity in 
most Earth orbit applications because it relates the rotating coordinate 
system, m which most measurements are made or related, to the inertial 
coordinate system, in which most calculations are made. This angle is 


( 1 - 8 ) 

The height 
values into 
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variously referred to as the Greenwich hour angle of the vernal equinox 
or as the Greenwich sidereal time — the latter is used herein. 

An approximate method for estimating the sidereal time for 0 hr 
GMT of any date is given as follows (more precise methods will be 
given following this short section). The sidereal time at 0 hr GMT is 
0 h on approximately September 22 (exactly zero on passage through 
the autumnal equinox) and increases at the rate of 3 m 56 s .5554 per 
day (see the next section “Time” following the coordinate discussion). 
Now, 3 m 56 s .5554 is very nearly equal to (4 m — 3 5 .5). Therefore, if N s 
is the number of days from September 22, then the sidereal time is 
approximately 


0g o % ( N s x 4) min - ( N s x 3.5) sec 

For example, find the approximate sidereal time at 0 hr GMT on 
March 1. September 22 is day 265, and March 1 is day 60. Then, 
the number of days from September 22 to March 1 is (365 - 265) 4- 
60 or 160 days. Thus, the sidereal time for 0 hr GMT on March 1 is 
approximately 

(160 x 4) min — (160 x 4) sec 

or 10 h 30 m 40 5 . The Astronomical Almanac for 1985 (AA85) gives 
10 /l 34 m 50 s . The reason for most of the difference is that in 1985, 
the autumnal equinox occurs at 6 hr GMT on September 21, which 
is 0.75 days earlier than assumed here. If we make this correction 
(which we usually wouldn’t do, as if we knew what correction to make, 
we wouldn’t have to use this approximate method), we would get 
10 h 33 m 37 s , which is just a bit more than a minute off. 

Incidentally, the day of the year for any given calendar date can be 
computed from the equations (Almanac for Computers 1980) 




(Nonleap year) 


(Leap year) 


> 


( 1 - 10 ) 


where M is the month (1-12) and D is the day (1-31). The symbol ( ) 
means “integer value of.” 

In order to compute the sidereal time to an accuracy of a few 
hundredths of a second, needed for many space applications, we need 
to introduce the concept of Julian date . 
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The Julian date associated with a given calendar date is equal to 
the number of days which have elapsed between noon, January 1 
4713 B.C. (November 24, 4713 B.C. in the modern calendar) and 
the date in question. This date was selected by astronomers as 
preceding the earliest recorded astronomical observations so that all 
known observations will have a positive Julian date, making it very easy 
to determine the time intervals between events. With all the sundry 
calendars which have been in use throughout history, the Julian date is 
generally found by referring to tables corrected to the present Gregorian 
calendar. However, for the years 1801-2099 A.D., a period probably 
inclusive for most of us, the conversion of Gregoric calendar date to 
Julian date can be carried out by the following formula (Almanac for 
Computers 1980): 


JD = 367F- 



UT 


+ 1721013.5 + — - 0.5 sgn (100F + M - 190002.5) + 0.5 


where 


( 1 - 11 ) 


Y = year (1801 < Y < 2099) 

M = month (1 < M < 12) 

D - day (1 < D < 31) 

and UT is the Universal time (Greenwich mean time, see later in 
this chapter for discussion). The symbol ( ) means “largest inte- 
ger function. The last two terms add up to zero for all dates af- 
ter February 28, 1900, so that these terms may be omitted for all 
dates. (There is another algorithm quoted by Blackadar 
(1984), which reportedly gives the correct Julian date for any calendar 

date from 4713 B.C. to 3500 A.D.-some readers might find this more 
useful.) 

Therefore, in order to compute the Greenwich sidereal time GST 
for any date and time T, use the following algorithm: 

_ J Se * UT = 0 in equation (1-11) and find the Julian date for 0 hr 

UT of the date in question. 
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Chapter 1 

2. Find the Greenwich sidereal time at 0 hr, 8 go , from the equa- 
tion (Escobal 1965, or the Explanatory Supplement to the Astronomical 
Ephemeris 1961, for example) 

6g o = 99°. 6909833 + 36000°. 76892T U + 0°. 000387087^ (1-12) 

(mod(360)) 

where the units are in degrees, and 

JD - 2415020.0 

Tu “ 36525. 

is the number of Julian centuries from noon on January 1, 1900. 

3. Let T be the number of minutes of mean solar time (usual clock 
minutes) from 0 hr UT to the actual time in question. Then, during 
this time interval, the Earth has rotated 

A 8 g = 0°.25068447T (l-14a) 

and hence the GST is 

8g = 9 go + A 6g (M4b) 

and now the Earth-fixed (rotating) coordinates can be related to the 
inertial coordinate system. Note the use of equation (1-12) implies a 
mean-of-date coordinate system, as nutational terms are not included 
in this equation. 

Both the mean sidereal time (precession only) and the apparent 
sidereal time (nutation also included) are tabulated in various almanacs, 
ephemerides, for example. Sidereal time is generally expressed in 
hours-minutes-seconds, even though it is usually used in degree units. 
To convert from decimal degrees to decimal hours, simply divide by 
15 — the Earth rotates 360/24 or 15 deg/hr. To display the accuracy 
of equation (1-12), tables 1-1 and 1-2 show the mean sidereal time 
computed from equation (1-12) and the seconds only of the mean 
sidereal time taken from the AA85. 2 The first of each month is shown. 
The maximum error in table 1-1, for the year 1985, is about 0.07 
second when comparing the almanac data with those computed from 
equation (1-12). The last column in the table shows the mean sidereal 
time computed from a set of formulas published in AA85 and given as 
follows: 

Q go = 100°. 4606184 + 36000°. 770057^ + 0°. 0003879337^ 

— 0°.000000258T^ 3 (1-15) 

2 The Astronomical Almanac 1985. 
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Table 1-1. Comparison of Computed and Published Mean Sidereal Times 
for First Day of Each Month for Year 1985 




Date 

Julian 

date 



Mean sidereal time, $ gt 

from — 

j Equation (1-12) 

AA85 

Equation (1-15) 



h 

m 

s 

S 

8 

Jan. 1 

2446066.5 

6 

42 

21.8975 

21.9674 

21.9674 

Feb. 1 

2446097.5 

8 

44 

35.1139 

35.1838 

35.1838 

Mar. 1 

2446125.5 

10 

34 

58.6641 

58.7341 

58.7340 

Apr. 1 

2446156.5 

12 

37 

11.8804' 

11.9505 

11.9504 

May 1 

2446186.5 

14 

35 

28.5413 

28.6115 

28.6114 

June 1 

2446217.5 

16 

37 

41.7576 

41.8279 

41.8278 

July 1 

2446247.5 

18 

35 

58.4186 

58.4889 

58.4888 

Aug. 1 

2446278.5 

20 

38 

11.6349 

11.7053 

11.7052 

Sept. 1 

2446309.5 

22 

40 

24.8512 

24.9216 

24.9216 

Oct. 1 

2446339.5 

0 

38 

41.5121 

41.5827 

41.5826 

Nov. 1 

2446370.5 

2 

40 

54.7284 

54.7990 

54.7990 

Dec. 1 

2446400.5 

4 

39 

11.3894 

11.4601 

11.4600 


r 

\ 


\ 


\ 


where 


T , _ JD - 2451545.0 
u 36525. 


which is the number of Julian centuries from January 15, 2000 A.D. 
Neither equation (1-12) nor equation (1-15) is exact. Equation (1-15) 
gives somewhat better results when applied to post-1984 data. The 
reason for this is that in 1984 the International Astronomical Union 
adopted a new revised set of physical constants which slightly changed 
some of the numerical constants in the time series — compare the 
numerical coefficients of T u and T' u in equations (1-12) and (1-15), for 
example. These new constants include the effects of general relativity 
and are arguably more precise than the older, pre-1984 constants, 
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Chapter 1 

although the absolute differences between the numerical results are of 
the order of hundredths of arc-seconds. 

Table 1-2, taken directly from AA85, shows both the mean (pre- 
cession only) and apparent (nutation included) sidereal times, and the 
differences between them, for the same dates. 

Table 1-2. Published Values of Apparent and Mean Sidereal Times for 
First Day of Each Month for Year 1985 



Apparent sidereal time 

Mean sidereal time, 
s 

Difference, 

9 

Date 

h 

m 

s 

Jan. 1 

6 

42 

21.1326 

21.9674 

-0.8348 

Feb. 1 

8 

44 

34.4222 

35.1838 

-0.7616 

Mar. 1 

10 

34 

57.9653 

58.7341 

-0.7688 

Apr. 1 

12 

37 

11.1409 

11.9505 

-0.8096 

May 1 

14 

35 

27.7665 

28.6115 

-0.8450 

June 1 

16 

37 

41.0078 

41.8279 

-0.8201 

July 1 

18 

35 

57.7705 

58.4889 

-0.7184 

Aug. 1 

20 

38 

11.0803 

11.7053 

-0.6250 

Sept. 1 

22 

40 

24.2781 

24.9216 

-0.6435 

Oct. 1 

0 

38 

40.8663 

41.5827 

-0.7164 

Nov. 1 

2 

40 

54.0502 

54.7990 

-0.7488 

Dec. 1 

4 

39 

10.7768 

11.4601 

-0.6833 


As seen, the differences are of the order of 1 sec, the mean sidereal 
time being about 1 sec later than the apparent sidereal time for this 
particular time interval. If the additional accuracy is needed (i.e., if the 
apparent sidereal time is needed), the following correction, called the 
Equation of the Equinoxes , can be added to the mean sidereal time as 
computed either from equation (1-12) or the sequence following it. The 
Equation of the Equinoxes is defined as the right ascension of the mean 
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Compilation of Methods in Orbital Mechanics and Solar Geometry 

equinox referred to the true equator and equinox. The expression for 
the apparent sidereal time is 

(d g ) apparent = ( 9 g ) mean + Aip cos e (1-16) 

where Atp, called the nutation in longitude, is given approximately by 
(Smart 1977, or Escobal 1968) 

Axb = - (17".233 + <y'.017T u ) sin 17 + 0".209 sin 2Q - l".273 sin 2 L 
— 0".204 sin 2 C ( 1 - 17 ) 

where Aip is in arc-seconds in this correction formula in which (Escobal 
1968) v 


fi longitude of mean ascending node of lunar orbit, measured in 
ecliptic plane from mean equinox of date, deg 
= 259°. 132750 - 1934° . 14200837^ + 0°. 002077787^ 

+ 0°. 00000222227^ q.j; 


L mean longitude of Sun, measured in ecliptic plane from mean 
equinox of date, deg 

= 279°. 6966778 + 36000°. 76892507’ u + 0° .0003025007^ (1-19) 


C — geocentric mean longitude of Moon, measured in ecliptic plane 
from mean equinox of date to mean ascending node of lunar 
orbit, and then along orbit, deg 
= 270°.4341639 + 481267°.88314177’ u - 0° .001 133337^ 

+ 0°.00000188897’3 m 2 0 


e = mean obliquity of ecliptic, deg 
= 23°. 4522944 - 0°.0130125T U - 0°. 00000163897^ 
-I- 0° .000000502787^ 


( 1 - 21 ) 


The correction to mean sidereal time using equations (1-16) to (1-21) 
is tabulated in table 1-3. 

The nutation corrections computed from equations (1-16) to (1-21) 
agree very well with the differences presented in the last column of 
table 1-2, the maximum error being about 0.014 sec in November. 
Escobal (1968, pp. 252—260 or pp. 304—305) gives another expression 
for the nutation in longitude which is reportedly accurate to 0.0001 
arc-sec, if such accuracy is needed. 
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Table 1-3. Computed Differences Between Mean and Apparent Sidereal Times 
for First Day of Each Month for Year 1985, With Intermediate Angular 

Values Included 


Date 

D, deg 

L, deg 

C, deg 

e, deg 

At p cos e, 
arc-sec 
(corrected) 

Jan. 1 

55.1508 

280.5969 

31.4281 

23.4412 

-0.8367 

Feb. 1 

53.5093 

311.1519 

79.8964 

23.4412 

-0.7360 

Mar. 1 

52.0266 

338.7501 

88.8355 

23.4412 

-0.7671 

Apr. 1 

50.3850 

9.3052 

137.3038 

23.4412 

-0.8125 

May 1 

48.7964 

38.8746 

172.5956 

23.4412 

-0.8539 

June 1 

47.1548 

69.4296 

221.0639 

23.4412 

-0.8243 

July 1 

45.5662 

98.9991 

256.3558 

23.4412 

-0.7222 

Aug. 1 

43.9246 

129.5541 

304.8241 

23.4412 

-0.6309 

Sept. 1 

42.2831 

160.1092 

353.2964 

23.4412 

-0.6443 

Oct. 1 

40.6944 

189.6786 

28.5843 

23.4412 

-0.7115 

Nov. 1 

39.0529 

220.2337 

77.0526 

23.4412 

-0.7344 

Dec. 1 

37.4643 

249.8031 

112.3445 

23.4412 

-0.6710 


The Greenwich sidereal time was defined as the hour angle of the 
vernal equinox 9 g (fig. 1-4). The local sidereal time is defined as the hour 
angle of the vernal equinox measured from the observer’s meridian H. 
Thus, from figure 1-4 it is seen that 


H = 0 g + \e = LST 


( 1 - 22 ) 


and, when expressed in tune measure, is the number of sidereal hours 
since the observer’s meridian was on the vernal equinox. Note that this 
quantity, 9 g + X E , was used in equation (1-2) which related r E in the 
rotating system to the vector r in the inertial axis system. 
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Compilation of Methods in Orbital Mechanics and Solar Geometry 

As stated earlier, the Greenwich sidereal time 6 g allows us to relate 
directly quantities calculated in the “inertial coordinate” system with 
those m the rotating, Earth- fixed axis system (fig. 1-4). 

In figure 1-6 we show the position of the vernal equinox and the 
Greenwich meridian, the two “z-axes” of coordinates, as well as the 
position of an Earth-fixed observer O and the spatial location of an 
Earth satellite. It must be emphasized that this picture is “frozen in 
time, as the Earth is rotating about the pole P, and the satellite is 
moving. However, at this instant, the satellite has a definite set of 
coordinates right ascension a s and declination <5 S , as measured in the 
inertial coordinate system, and a definite latitude ip s = 6 S and longitude 
A s measured in the Earth-fixed coordinates. The subsatellite point S is 
defined to be the point on the surface when the geocentric radius vector 
of the spacecraft (S/C), r s , pierces the Earth’s surface. The angular 
coordinates of S are those of the spacecraft. 

Let ft* be a vec tor to the observer in the rotating coordinates 
lV>o,A 0 ); He is the magnitude of the Earth’s radius, 



Figure 1-6 Sketch relating position of observer in rotating Earth-fixed axes to 
inertial position of celestial body. 1 
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The coordinates of the spacecraft are similar, 


r s = r s 


cos xp s cos 
cos rp s sin A s 
sin t p s 


(1-24) 


These two vectors form a plane, as shown in figure 1-7. The angle 
z is called the zenith distance, where zj is the zenith distance of 
the spacecraft as measured at the center of the Earth, and z 0 is that 
measured by the observer O. For near-Earth objects, such as an Earth 
satellite, the Moon, and for some purposes the Sun, we have zj- ^ z Q 
(parallax effect). For observation of stars, though, as r s approaches oo, 
zj* approaches z 0 , and z 0 can easily be found from 


_ Re ' r s 
cos 20 “ 

= sin xp Q sin xj) s + cos xp D cos x p s cos (A s — A 0 ) (1-25) 


Define the radius from the observer to the spacecraft, 

p = r s - Re 


and hence 


r 


COS Zj- — 


Re P 

ReP 
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Compilation of Methods m Orbital Mechanics and Solar Geometry 
which can be put into the form 


cos 


z T = 


r s cos Zq - Rg 

foe + r i ~ 2Rgr s cos z c )^ 2 


(1-26) 


From equation (1-26), it can be seen that, as r s approaches oo, z„ 
approaches z x as expected. 

The sine and cosine of the azimuth of the spacecraft, A, can be 
found by applying the law of sine and the law of cosine for sides to the 
spherical triangle POB of figure 1-6 


sin A = — ^ sin ~ A °) ' 

sin z 0 

cos A = — ** - sin cos Zo I ^* 27 ) 

cos ip 0 sin z 0 

Don t forget a s , 6 S , z 0 , zj - , and A are time-dependent quantities. 

Another coordinate system that is very often encountered is a local 
axis set fixed in some definable way to the orbiting spacecraft. Mea- 
surements with on-board instruments and/or measurement direction 
vectors are made in this local system, and one must frequently trans- 
form vectoral quantities back and forth between this local system and 
either the inertial system or the Earth-fixed system defined earlier. 

The most fundamental local system is one which is defined com- 
pletely in terms of the dynamic variables r and v, the position and 
velocity vectors of the spacecraft given at some time t in the inertial 
coordinate system. These can perhaps best be visualized by thinking 
of the spacecraft as an airplane flying in the “normal” flight position 
The unit vector 


e r = 


r 

r 


ti 

f-2 

^3 


(1-28) 


defines a unit vector in the “up” or local vertical direction, positive 
when pointing away from the center of the Earth. 


A second vector 


v x r 

\v x r | 


mi 

m 2 

m 3 


(1-29) 


is a unit vector pointing out the “right wing” of the spacecraft and 
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Chapter 1 

is opposite in direction to the orbital angular momentum vector. The 
third vector 


e v — c r x e n — 


n l 

n 2 

n 3 


(1-30) 


is a unit vector pointing “forward.” If the orbit were exactly circular, 
then e v would be parallel to the velocity vector v. 

These three mutually orthogonal unit vectors are sketched in 
figure 1-8. 



Figure 1-8. Sketch showing relation of spacecraft-fixed unit vectors e r , e n> and 
e v to inertial axes. 


or in matrix form 


= 


are simpl) 


V n 


Vv 


V r 

v n 

r 

V v 

— 

m Vr m 



h h ^3 
m\ m 2 m3 
ni ri2 ri3 



Vlx 


v Iy 


yiz. 


(1-31) 


(1-32) 


The matrix of equation (1-32) is a pure rotation matrix, and hence 
transformation from the local S/C axis system to the inertial axes is 
readily accomplished by using the matrix transpose. In particular, let 
v L be a vector defined in the local S/C system by its magnitude and by 
the azimuth angle Ai and elevation angle n as shown in figure 1-9. 
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Compilation of Methods m OTbital Mechanics and Solar Geometry 

The local vector components are 

v n — vi cos 7 1 sin A l 1 

V v = v L cos 7 L cos A l i (1-33) 

v T = v L sin 7 L J 

and the components of this vector in the inertial system are then 


(1-34) 


v Lx ^1 

mi 


v Ly = h 

m2 

n 2 

v Lz\ ^3 

m 3 

n 3 

e r 





Figure 1-9. Description of orbiting vector v L in local spacecraft axes defining 
azimuth A i and elevation angle & 

Equation (1-2) then provides the link between the inertial coordi- 
nate system and the Earth-fixed coordinate system, and hence, vector 
quantities can easily be transformed between the Earth-centered coor- 
dinate and the local S/C system. 

Other instrument or spacecraft-specific coordinate transformations 
can, of course, be readily defined relative to the local dynamic coor- 
dinate system defined by equations (1-33). One has merely to chain 
the sundry transformation matrices together to provide transformation 
links between any of them and the inertial or Earth-fixed systems. 

Time 

The concept of time is one that most of us take pretty much for 
granted. But time is one of the most difficult concepts to resolve in 
orbital mechanics and astronomy because of the many very small effects 
that creep into its measurement and definition. As precision methods 
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were developed, many time concepts have had to be redefined because 
parameters which were once thought to be constants or absolutes were 
found to be indeed variables. Entire books have been written on 
this subject alone. Consequently, the following discussion will treat 
only the main features of time and will serve only to permit the 
gross features of and differences between sidereal time (or time with 
respect to the fixed stars) and solar time (time with respect to the 
Sun) to be understood and appreciated. Chapter 10 of Green (1985) 
contains a fairly thorough discussion of the various times considered in 
modern astronomy. However, these concepts are generally much more 
stringent than those required in most Earth orbital applications of orbit 
mechanics theory. (See also Smart 1977, chap. VI.) 

In figure 1-10 is shown a portion of the Earth's orbit around the 
Sun. Suppose that, at point (T), there are two observers on the Earth 
located exactly 180° apart in longitude. Observer A is watching the 
stars and observer B is watching the Sun. Suppose further that they 
are in constant communication with each other (relativity effects are 
ignored) and that at the exact instant that observer A reports a star on 
his meridian, observer B reports that the center of the Sun is exactly on 
his meridian. Now, observer A has been watching the stars for years. 

He has constructed a clock, based on stellar time, which is extremely y- 

(infinitely) accurate. He has defined 24 hr of sidereal time as the interval 

between two successive passages of the same star over his meridian (see 

the note at the end of this chapter). Each hour is divided into 60 

sidereal min and each minute into 60 sidereal sec. Observer A has set 

and calibrated an identical clock which he has given to B. So, at the 

instant ®, both A and B start their clocks. 



Figure 1-10. Position of Earth relative to Sun on two consecutive days, ® and 
and 6 months later, 
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At point © (fig. 1-10), which occurs one “day” later, observer A 
calls out that his star is again on his meridian and that his clock reads 
, T* - ° b ® erver B notices that, sure enough, 24 hr has passed on his 
clock too, but the Sun has not yet appeared on his meridian. In fact, he 
finds that he must wait an additional 3 min and 56 sec for the Sun to 
fine up m his instrument. The following day, he finds that he must wait 
7 mm 52 sec, etc., with each day being about 3 min 56 sec longer than 
the day before. At the end of one full revolution of the Earth around 
the Sun, 1 yr, observer B finds that his observation of the Sun’s center 
on the meridian is exactly 24 sidereal hr late. 


When the Earth arrives again at point ® , 1 yr later, observer B 
would find that he has made 365.2422 revolutions with respect to the 
Sun. Observer A, of course, would find that he has made 366 2422 
revolutions with respect to the fixed stars. (When B makes one daily 
rotation with respect to the Sun, A has made one full revolution with 
respect to the stars plus a little more — in fact, 3 m 56 s more or a in 

Jj*?® I'™, TheSC “ littIe mores ” add U P t0 exactly 1 full day in 1 yr 
If the Earth were not rotating, the Sun would still make one apparent 
revolution around the Earth in 1 yr. 


Now, observer B finds that the amount he must wait each day for 
the Sun to appear on his meridian is not exactly 3 m 56 s every day. The 
interval is somewhat longer in December and shortest in September (see 
hg. 1-11), but the mean, averaged over 1 year, is 3 min 56 sec 


1440 min/day 
365.2422 


3 m 56®.5554 


Observer B correctly attributes these daily variations to two factora- 
ge apparent orbit of the Sun around the Earth is not a circle but an 
ellipse and the ecliptic, the plane in which the Sun appears to move, is 
mclined to the equator. Both these factors cause a nonunifonn motion 
of the Sun about the Earth (fig. 1-11). 

«ftj^ gUre 1 ' 11 . shows a P lot of the difference between the solar day and 
86400 mean solar sec for each day of the year. The solar day is longest 
in late December for two reasons: (1) the Earth is near perihelion and 
has the largest angular velocity in its orbit and (2) the Sun is also near 
the winter solstice, and consequently, is moving essentially parallel to 
the equator. The solar day is shortest in mid-September (and almost 
as short in mid-March) because the Sun is crossing the equator at these 
times and hence has the smallest component of its angular velocity 
projected onto the equator. 
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Month 


Figure 1-11. Plot of quantity At = Length of solar day — 86400, in mean solar 
seconds, versus day of year. 


Observer A now asks the very legitimate question: Why can’t B use 
A’s clock to keep his (B’s) time? The answer is, of course, that he can, 
but B would not find this aesthetically pleasing and it would, in fact, 
cause him considerable “headaches” later on. If B used A’s clock, then 
at the starting position shown at © in figure 1-10, B would define his 
time as local noon and A’s time as local midnight. A new day would 
start for A and would start for B 12 hr later. However, 6 months later, 
A and B would be at position (3) of figure 1-10. Now, keeping A’s time, 
midnight for B would occur when the Sun is directly overhead on his 
meridian, that is, at midday. Also, during the 6 months which have 
elapsed, if B retains the notion to start a new day 12 hr after the Sun is 
directly overhead for him, he would find that the change from one day 
to another would occur at various times of day throughout the year. 
This would certainly lead to some confusion if not some real practical 
difficulties. 

Long before the beginnings of recorded history, man has regulated 
his everyday affairs by the Sun — working, hunting, etc., by day and 
sleeping at night. Since even now, most people work during the day, 
it is convenient to have days change from one to the other during the 
hours of darkness, that is, at midnight. Therefore, a time geared to the 
Sun would be extremely useful. However, the direct use of the observed 
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Sun would produce days which vary in length throughout the year for 
the reasons mentioned earlier. 

In order to construct a clock which keeps uniform solar time of 
some sort, B proposes the use of a fictitious Sun, the fictitious mean 
Sun, which orbits the Earth in the equatorial plane, rather than the 
ecliptic, and moves at a uniform rate such that the fictitious mean Sun 
completes one revolution about the Earth in exactly the same time 
period as does the true Sun i.e., in 365.2422 mean solar days. Since 
the fictitious Sun moves in the Earth’s equatorial plane at a uniform 
rate, its right ascension is increasing at a uniform rate. The difference 
between the right ascensions of the mean Sun and the true Sun is called 
the equation of time (see fig. 1-12): 


E? = RAMS - RATS (1-35) 

and can readily be computed from simple orbit mechanics (see, for 
example, Smart 1977) and, to terms of the second order in the Earth’s 
eccentricity, is given in units of radians by 

Et = V sin 2 L - 2e c sin M s + 4t/e c sin M s cos 2 L 
1 O r 5 O 

- -y sin 4 L - -ei sin 2 M, + . . . (l- 36 ) 

where 


y = tan 2 


£ 

2 


(1-37) 


£ = obliquity of ecliptic, equation (1-21) 
e e = eccentricity of Earth orbit 

= 0.0167514 - 0.00004187^ - 0.0000001267^ (1-38) 

and where the right ascension of the mean Sun, RAMS, and the Sun’s 
apparent mean anomaly (see chap. 4 for the definition of some of the 
orbital element concepts given here) are given by (Escobal 1968) 


L — RAMS (eq. 1-19) or L + A if cos e if correction 
needed 

M s = 358°. 475844 + 35999°.049757 u - 0°.000157’2 
- 0°. 000003337^ 


(l-39a) 

(l-39b) 
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Figure 1-12. Sketch showing geometrical relationship between true Sun and mean 

Sun. 

The nonlinear terms in this equation are due to the easterly motion 
of the Earth’s perihelion. 

See chapter 10 of Green (1985) for a modern definition of Ej- and 
an equation for its calculation. The concept of the equation of time is 
no longer used in modern astronomy, although as a computational tool, 
it is still quite useful. The fall from grace of this concept came about 
when it was discovered that the rotational speed of the Earth about its 
axis was not truly a constant. The time definitions we have given are, 
therefore, only approximately correct, as they hinge on the constancy of 
this quantity. Equation (1-35) is still correct, however, if ephemeris time 
is used. Ephemeris ’time is the time used in the differential equations 
of motion. It is thus based on a dynamical concept, rather than the 
geometrical concepts used up to now. 

When B sets his clock to the fictitious mean Sun, he is now 
measuring mean solar time, or more commonly, the ordinary clock time, 
with which we’re all familiar. The international time, universal time, 
formerly called Greenwich mean time, is computed a s mean solar time, 
with Greenwich mean noon being the instant that the mean Sun is on 
the Greenwich meridian. 

As pointed out, many of the simplified time concepts discussed here 
are only approximately true because they depend on the constancy of 
the rotational rate of the Earth with respect to the fixed stars. Within 
the last half-century or so, it was found that this premise is not strictly 
true, and other definitions of time, which do not depend on the Earth’s 
rotation, have been introduced. These very small differences, however, 
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Compilation of Methods in Orbital Mechanics and Solar Geometry 

are generally of little importance in applications of these equations to 
Earth orbit mission design and studies and are generally ignored. These 
concepts are quite accurate enough for determining the position of the 
Sun, for example, as will be shown in chapter 5. Most of the new texts 
on Spherical Astronomy address this question of time most adequately 
(See, for example, Green 1985, Taff 1981, Taff 1985, and the edition of 
Smart s classical text revised by Green listed herein as Smart 1977.) 

In practice, of course, it would be impractical for every location 
on the surface to keep its own mean solar time. To standardize this 
somewhat, the Earth was, by international agreement, divided into 24 

tin !fo Z T, eS ° f a PP roximatel y 15 ° of longitude. The time zone centered 
at 0 of longitude, the Greenwich meridian, and extending for 7.5° on 
either side^of the 0° meridian keeps Greenwich mean time. The time 
zone at 15°± 7.5° west is the first time zone, and so on. For the United 
States, eastern standard time is referenced to the 75th meridian, central 
time to 90°, mountain time to 105°, and Pacific standard time to 120° 
west longitude. Thus, since 75° W corresponds to 5 hr of time we can 
convert standard times in the four time zones to GMT as follows: 


r Eastern standard time 

+ 

5 hours ' 

Central standard time 

+ 

6 hours 

Mountain standard time 

+ 

7 hours 

» Pacific standard time 

+ 

8 hours ; 


= GMT 


(subtract 1 hour from these numbers for daylight savings time). 

With all the perturbations and undulations which time can take, th< 
really important thing to remember for our purpose is that a siderea 

366.2422 mean sidereal days, and the mean solar year has 
odo. 2422 mean solar days. Thus, 


r 


24 hours of mean solar time = 24 x j^ 6 ' 2422 

365.2422 


24 3 56 a .555 of mean sidereal time 


and 


24 hours of mean sidereal time = 24 x — 5,2422 

366.2422 


— 23 56 4 s . 091 of mean solar time 


The ratio 365.2422/366.2422 is, of course, the ratio of any (mean solar 
time interval)/ (mean sidereal time interval). This is the source for 
the weird-looking constant 0.25068447 in equation (1-14). The Earth 
rotates 0.25° in 1 sidereal min. Therefore, in 1 mean solar min it rotates 
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366.2422 , , , 

0.25 x = 0.25068447 deg/mean solar min 

365.2422 

Note: the scenario between the two astronomers was , of course, 
highly simplified in many ways, and in fact , one concept used was 
intentionally erroneous at that point for clarity . A sidereal day is 
actually defined as the time interval between two successive passages 
of the vernal equinox over the same meridian, rather than the interval 
between two successive passages of a fixed star ( one whose proper motion 
is essentially zero). Since, as mentioned earlier, the vernal equinox 
is moving westward at nominally 50 arc-sec /yr, the sidereal day is 
actually a bit shorter than it would be if the fixed star were used for 
time definition. The precessional constant is 50.2564 circ-sec/yr and 
is measured along the ecliptic. Its component along the equator is 
thus 50.2564 cos 23.44 = 461091 arc-sec/yr, or 0.12624 arc-sec/mean 
sidereal day, and hence, the day as defined is 0.12624/15 or 0.008416 
sidereal seconds shorter than it would be if it were defined relative to 
the “fixed” stars. This amounts to about 1/120 sec. (See, for example, 
Motz and Duveen 1966.) 
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Chapter 2 

Gravitational Field of the Earth 

The external field of the Earth can be described mathematically 
in many ways. Because of its rotation, and its plasticity, especially 
during its early formative years, the Earth is very nearly an oblate 
spheroid. Hence, for many analytical purposes, it is convenient to ex- 
pand the Earth’s external gravity field in a series of spherical harmonics 
(Heiskanen and Moritz 1967, p. 59): 


n=0 m=0 L 


Pnm (“07 


r n4- 1 


B n 


S n m (0j A) 1 


r n+ 1 


( 2 - 1 ) 


in which 0 and A are the latitude and longitude, respectively, and r is 
the distance from the origin to the point at which the potential is to be 
computed. The constants A nm and B nm are integrals determined by 
the internal mass distribution of the Earth, and 


B n m (07 A) = P nm (sin 0) cos mX (2-2) 

S nm (0, A) = Pnm{ sin 0) sin mA (2-3) 


in which P nm are associated Legendre polynomials, are the spherical 
harmonics. 

The integrals A nrn and B nrn are in practice determined from the 
precision tracking of large numbers of Earth satellites and then per- 
forming statistical fittings to inverted data to determine the values for 
these constants, which best fit the observed tracking data. The God- 
dard Space Flight Center (GSFC) has done much of this work to high 
precision. 

Equation (2-1) essentially assumes that the Earth is a distorted 
sphere. The n — 0 term is the spherical part of the gravity field (the 
inverse square part), and the subsequent terms are needed to describe 
the departure of the body from sphericity — the greater this departure, 
the more terms are needed to describe the external potential to a specific 
accuracy, and the larger are the constants A nrn and B nm . 
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The centra] force part of the gravity field is usually described by the 
term 


V cf = 7 (2-4) 

where /z is the gravitational constant. If one takes the n = 0 terms out 
of equation (2-1) and writes successively 


n=l m= 0 L J 

{ °° n n \ 

1 + (~) [^nm #nm(^, A) + /f nm S nm (^ } A)] l (2-5) 

n= 1 m=0 J 


where 


^OO = M 


P n T 

J nm 


R?K n 


^nm 

^oo 

Bnm 

-4 00 


and R<. is the “equatorial radius” of the body, then equation (2-5), as 
written, describes the gravity field as a central force term (the term fi/r) 
plus a “perturbation potential” (the summation terms of eq. (2-5)). 

If the origin of the coordinate system coincides with the center of 
the mass of the central body, a usual assumption, then (Heiskanen and 
Moritz 1967, p. 63) 


Jio = J\1 = K\1 = o 

and we write equation (2-5) as 

{ oo n n \ 

1 + (~) I ^ nm ^m(^, A) + K n m Snm(V>, A)J 1 (2-6) 

n=2m=0 j 

The various types of harmonics in equation (2-6) are sketched in 
figure 2-1 and are identified by 

m = 0, zonal harmonics 
m = n ^ 0, sectorial harmonics 
m ^ n ^ 0, tesseral harmonics 
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Zonal 

fnO 



Sectorial 
Ain 
(" *0) 


Tesseral 
F nm 

(n ^ //I / 0) 




Figure 2-1. Sketch showing various types of harmonic coefficients. White area 
implies elevation above; black, elevation below mean spherical surface. 


The zonal harmonics are functions of latitude only and hence can only 
reflect north-south variations in the gravity field. 

Many of the precision orbit prediction programs used by GSFC, 
NOAA, and others use large numbers of terms from equation (2-6). 
GSFC, for example, has (at least) two such programs, one of which goes 
up to n = m = 8 and the other to m = n = 21. The GEM-8 model is the 
one used for SAGE I and II and SAM orbit work (these programs include 
other than gravitational effects, e.g., drag, atmosphere models which are 
affected by sunspot activity, relativistic effects, other planetary effects). 

For many applications it is generally found that the sectorial and 
tesseral harmonic terms in equation (2-6) can be neglected and that only 
the first few zonals are needed. Thus, if we retain only the m = 0 terms, 
then all the K nm terms vanish, and letting J n o = Jn, equation (2-6) 
can be written as 


v = * 


r 



n 

J n P n (sin ip) 


(2-7) 


where P n is the nth-order Legendre polynomial. The first few of these 
are (Heiskanen and Moritz 1967, p. 23) 


P 0 (x) = 1 

Pl{x) = X 


P2(x) 
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Compilation of Methods in Orbital Mechanics and Solar Geometry 
and the rest can be generated from the recursion formula 

p„w = + e2±i>xj»„ (I) (2 . 8) 

n n \ \ / 

If we go as high as J 6 , we can write equation (2-7) as (Escobal 1965) 

V= r{ 1 + ^( : 7 i ) (1-3 sin 2 VO + y (y-J (3-5 sin 2 iP) 

x sin iP - y (3 - 30 sin 2 iP + 35 sin 4 VO ~ y (~ ) 5 

x (15 — 70 sin 2 ip + 63 sin 4 VO sin ip 


J 6 / J ^ \ D X 

+ i6 It) _ 105 sin2 ^ + 315 sin4 ^ - 231 sin6 vo + • • • > 

- (2-9) 

The constants in equation (2-9), as obtained from GSFC in March 1986, 


are 


H = 398600.64 km 3 /sec 2 
R = 6378.14 km 
J 2 = +1082.6271E-6 
J 3 = -2.5358868E-6 
J 4 = -1.6246180E-6 
J 5 = -0.22698599E-6 
Je = +0.54518572E-6 

The potential (fq. (2-9)) and these constants are used in some of the 
author’s orbit programs to generate short-term ephemeris data for use 
in the SAGE II and SAM II data reduction. (See chap. 4.) 
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Chapter 3 
Orbital Elements 

As will be seen in chapter 4, the equations of motion in Cartesian 
coordinates of a satellite orbiting about an oblate central body (zonal 
harmonics only) are rather simple to write, and with modern computers, 
it is possible to numerically integrate these rapidly and accurately. 

However practical for number generation, Cartesian coordinates are 
not the most useful coordinates for visualizing or otherwise describing 
a spacecraft orbit. A time sequence of spacecraft position and velocity 
vectors by itself has little pictorial value and, consequently, conveys 
little information about the evolution of the orbit. 

There are several sets of “orbital elements” used in astronomy, 
astrophysics, space sciences, etc., each of which is most useful in 
the specialized application for which it was conceived. It takes six 
independent coordinates to completely specify the state of an orbiting 
spacecraft and permit the determination of its future (or past) state (for 
example, three position and three velocity Cartesian components), and 
hence, it also takes six independent orbital elements. For Earth orbit 
analysis or “Keplermanship” (Escobal’s term for playing games with 
the two-body equations), the set described below is the one favored by 
this writer, in both its utilitarian and interpretive senses. 

From the first of Kepler’s laws, we know that (for central body 
motion) the orbit of an Earth satellite is an ellipse with the center of 
the Earth at one of the foci. The ellipse has two axes — the major axis 
AB and the minor axis CD in figure 3-1. The origin is at point O, the 
center of the Earth. The spacecraft S is located by the radial position 
r and the angle / measured from the major axis and where / = 0 is 
defined by the radius OA where the spacecraft is closest to the Earth, 
the perigee. The equation of this ellipse is 

r = , a{1 ~ e — f (3-1) 

1 + e cos / 
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Figure 3-1. Geometry of ellipse in plane. 


where the two orbital elements , 


a = semimajor axis, km 
e = eccentricity 


determine the size and shape of the ellipse. 

For central force motion, the orbit lies in a fixed plane, passin, 
through the center of the Earth. Two additional orbital elements 


H — right ascension of ascending node, deg 
i = inclination to equator, deg 


locate the position and orientation of this plane in the “inertiaT 
coordinate system defined earlier. (See fig. 3-2.) Note that by genera 
agreement, if * < 90°, the orbit is called a “posigrade” orbit, whereas 
for 90 <i< 180 , the orbit is referred to as “retrograde.” 

The fifth orbital element, 


w — argument of perigee, deg 

(fig. 3-2), locates the position of the major axis of the orbit in this 
plane (more specifically, the position of perigee) and is measured in the 
direction of motion from the ascending node of the orbit. 
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Chapter S 



Figure 3-2. Sketch defining location of orbital plane and perigee point in inertial 
coordinate system. 


The sixth element must somehow relate the position of the space- 
craft in its orbit to some specific time. The angle /, called the true 
anomaly , is awkward to use for this purpose. According to the second 
of Kepler’s laws, the radius vector of the spacecraft sweeps out equal 
areas in equal time increments. Since the spacecraft is closer to the 
Earth at perigee (/ — 0) than it is at apogee (/ = 180), the spacecraft 
must move faster at perigee than it does at apogee. This means that the 
angular rate df jdt is not constant around the orbit. A mean angular 
rate of the spacecraft can be shown to be 


2n nr 

Period V a 3 


(3-2) 


and is a cons taint of the orbit for a central force. Let t Q be the time the 
spacecraft last passed through perigee. Then the angle 


M = n(t — t Q ) (3-3) 

is an angle which hats the desirable property of increasing at a uniform 
rate. This amgle is called the mean anomaly amd is selected here as 
the sixth orbital element. In order to relate the mean anomaly to the 
true anomaly, we must first introduce a third amomaly, the eccentric 
anomaly E , which relates / to M as follows: 
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, E /T—T / 

tan 2 = v TTi tan 2 M 

and Kepler’s equation 

M = E - e sin E (3-5) 

These transcendental equations are still not directly solvable (i.e., to 
find / given t). However, for orbits with small eccentricity (near-circular 
orbits) as most near-Earth orbits are, the following series (Smart 1977) 
is useful — angles are in radians: 


( e 3 \ e 2 

E = M + [e - — sin M + — sin 2 M 

\ C 2 

3 e 3 

+ — sin 3 M + . . . 

O 


/ = 2e- -J 

+ sin 3 M 

12 


be 2 

sin M -i — — sin 2 M 
4 


(3-6) 


(3-7) 


If the eccentricity is large, or if extreme accuracy is needed, the Taylor 
expansion of equation (3-5) gives the iteration formula 


En + 1 + En ~ 


Mj' — M n 


1 — e cos E n 

in which, given the true value of M T , find E n from equation (3-6) and 
M n from equation (3-5). Then equation (3-8) gives a corrected value 
E n+ 1- For a second iteration, set E n = E n+X from the last iteration, 
and repeat to convergence. This procedure generally converges to six 
or seven decimals in three or four iterations, even for values of e near 
unity. 


There are two basic advantages for using orbital elements. First, the 
orbit is much easier to visualize, as the orbital elements describe the 
total geometry of the orbit — its size, shape, and orientation in space. 
Second, and perhaps more important, if the Earth were a perfect sphere 
(i.e., the gravity field describable by the inverse square law only) and 
there were no other external forces acting on the spacecraft, then five of 
the six orbital elements o, e . t, fl, and uj — would be constants , only 
M would change with time, and this change would be a simple linear 
one. Even for a spherical Earth, all six of the Cartesian components 
change continually and dramatically with time at each integration step. 
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It seems reasonable to expect, then, that if the gravitational poten- 
tial field of the central body only deviates very slightly from that of a 
sphere as does the Earth, the orbital elements a, e, i, ft, and u) would 
change very slowly with time, if at all, and these changes would be 
predictable to varying orders analytically. This is, of course, precisely 
what happens for an Earth satellite. Very detailed and complex analyt- 
ical studies (e.g., Kozai 1959, Brouwer 1959, King-Hele 1958, Garfinkel 
1959, Merson 1961) show that it is possible to account for practically all 
the oblateness perturbative effects of the Earth with a single term (^ 2)1 
and to first order it is found that a, e, and i are constants, and ft and u 
change linearly and slowly with time (order of a few degrees per day, or 
less, depending on the inclination of the orbit). These changes are easy 
to compute (Escobal 1965 or McCuskey 1963). After the equations of 
motion are introduced in chapter 4, a brief series of numerical exam- 
ples will show the effects of various harmonics on both the Cartesian 
components and on the orbital elements. 

It is, of course, possible, if not imperative, that we be able to transfer 
from Cartesian coordinates to orbital elements and vice versa. The 
following algorithms, applicable to circular and/or elliptical orbits only, 
are used for this purpose (see, e.g., Escobal 1965 or McCuskey 1963). 

1. Orbital elements to Cartesian coordinates (CONCAR) 

(a) Find E from (M,e) as described earlier (eq. (3-6)) 

(b) Then 


x u = a (cos E — e) 

Zu = a VT — e 2 sin E 

Xu — -aE sin E 

z u = aE y/l- e 2 cos E 

where 

n 

E = 

1 — e cos E 

(c) Compute the unit vectors 


P = 


cos u cos ft — sin lj sin ft cos i 
cos uj sin ft + sin a ; cos ft cos i 
sin a; sin i 


T 


Q = 


— sin uj cos ft — cos u sin ft cos i 

— sin uj sin ft + cos a; cos ft cos i 

cos u) sin i 
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(d) Then compute the spacecraft position and velocity vectors 
from 

r = XuP + y^Q 
r = XojP + y^Q 

2. Cartesian coordinates to orbital elements (CARCON) 

(a) From angular momentum vector 


yz- zy 


h sin i sin f2 

zx — xz 

= 

—h sin i cos fi 

xy -yx _ 


h cos i 


Find h, t, and fi 

(b) With r — |r|, v = |r|, find the semimajor axis a from 

1 _ 2 _ v 2 
a r fj, 

where for the Earth, fi = 398600.64 km 3 /sec 2 . 

(c) Find the eccentricity 



(d) Define u = w + /, then 


r cos u = x cos fi -f- y sin Q 
2 

r sin u = 

sin i 

which gives the proper quadrant for u 

(e) Find the true anomaly / 

e cos f = - — 1 
r 

r V l r ' * 
sm / = — | — _ J 

where p = a( 1 - e 2 ) 

(f) Finally, then, 


r ^ 


40 


u> = u — f 

(g) Find E from equation (3-4) and M from equation (3-5). 
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This completes the orbital element set. 

Note that these equations even hold for circular orbits (e = 0), if we 
agree to set oj = 0 and measure / (or M, E ) from the ascending node 
of the orbit. 

The algorithm, CONCAR, changing from orbital elements to Carte- 
sian coordinates, is quite useful in determining time histories of the 
Cartesian coordinates, which parameters are frequently used in deter- 
mining other geometric parameters and properties in which the mission 
design engineers might be interested. Given the set of orbital elements 
at time t 0 , the mean anomaly M is found at t + At from equation (3-3) 
as 


M = M 0 + n At 

The eccentric anomaly is found from equation (3-6), as described earlier, 
an d then the Cartesian state vector is found from the algorithmic 
equations. This process can be repeated at any time interval A t, and 
the results hold as long as only a central force field (all J n — 0) can be 
assumed. In actual practice, this time might be of the order of 10 min 
or so, and then the small perturbations begin to make themselves felt, 
and another algorithm, to be described in chapter 4, must be used. 
This requires using the analytical expressions for u and fi explained in 
chapter 4 to continuously update u> and Q, thus enabling the algorithm 
CONCAR to be used over extended time periods. 

A third algorithm, especially useful if one wants to preserve the 
identity of the initial Cartesian coordinates, is the so-called /- and 
^-series method. (See, for example, Escobal 1965.) Here, the radius 
vector r{t) and velocity vector ir(t) at any time t are expressed in terms 
of the initial values of the position and velocity vectors, and the /(£, t 0 ) 
and g{t,t Q ) parameters 

f(t) = /(*, to) r Q + g(t, to)r 0 (3-8) 

r{t) = /(Mo) r o + g{t,to)r 0 (3-9) 

where r Q and r 0 are the initial position and velocity vectors, respectively. 
The / and g terms were originally derived as an infinite series involving 
the zeroth, first, and second derivatives of the position vector as 
const ant coefficients and the time as the independent variable. However, 
for central force motion, these series can be expressed in analytical form 
as 


f(t, to) = 1 - - [1 - cos (E - Eo)} (3-10) 

r Q 

41 





t f . * i 


\ 7 


_V ' V V - ir 1/ l .'. I ’■ * ' 


V J * ’ -17- 


l A 





Compilation of Methods m Orbital Mechanics and Solar Geometry 

g(t, t a ) = t- [( E - E 0 ) - sin ( E - E 0 )\ (3-11) 

and 


f (t, to) — - 

g{t,t 0 ) = 1 

in which 


r = a(l - e cos E) (3-14) 

The parameters a, e, and E a are computed from algorithm 2 
(CARCON) at t = 0, and E at any subsequent (or previous) time 
can be found from equation (3-6) with the help of the iteration process 
described previously. This procedure offers no real computational ad- 
vantage over that of algorithm 2 alone. It has, in fact, the disadvantage 
of not readily incorporating the oblateness effects through the fi and u 
terms as described earlier. It does, however, retain the identity of the 
initial Cartesian state, as mentioned, and is quite useful in studying 
the effects of errors in the initial conditions, as the variance-covariance 
matrix of errors in the values of future coordinates can be constructed 
directly from equations (3-8) and (3-9). 

As pointed out earlier, five of the six orbital elements would be 
constant if the Earth were a perfect sphere and there were no other 
outside perturbations acting on the orbit. The effects of the nonspher- 
lcal components (in the spherical harmonic sense) of the gravity field 
the perturbation forces, are very small for a typical Earth satellite, and 
hence, the orbital elements, instead of being constant, vary very slowly 
in time. The most elaborate analyses referenced earlier (Kozai 1959 
and the others) indicate that these changes are combinations of the 
following: 

Secular changes: these are linear, or at best quadratic, changes, which 
always proceed in the same direction. The elements Q, u>, and M 
are the only ones which experience secular changes. 

Long-period terms: periodic changes of small amplitude whose periods 
are of the order of 80—100 days and longer. 

Short-period terms: periodic changes of small amplitude whose periods 
are small (integer or half-integer) multiples of the orbital period. 

These terms are sketched in figure 3-3. 
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Escobal (1965, chap. 10) discusses perturbations of the orbital 
elements thoroughly, as does the classical text by Moulton (1914). A 
typical orbital element can be most generally represented by an equation 
of the form 

q — q 0 + q{t - t 0 ) + K\ cos 2cj 4- K .2 sin (2/ + 2ui) + H.O.T. 
where 

q ~ secular term 
q Q = mean value of the element 
K\ = long-period term 
if 2 = short-period term 

Higher order terms (H.O.T.) in other multiples of w and / are, of course, 
also present, and any of the constants can be zero. 

In any event, as a result of these perturbations, the orbit of an 
Earth satellite is never a true ellipse. However, as seen in the preceding 
algorithms, the specification of the spacecraft state vector permits the 
computation of a unique set of orbital elements, and thus at each instant 
of time, a unique set of orbital elements can be associated with the orbit. 
If at a given instant of time all the gravity perturbations were turned 
off and only the central force part of the gravity field were allowed to 
remain, the spacecraft would then truly be in the elliptical orbit defined 
by the orbital elements at that moment. This constantly changing set 
of orbital elements is referred to as the set of osculating elements and is 
specified at a unique time. Point a at time t in figure 3-3 identifies the 
osculating element at this time. It has been found through many sets 
of calculations that, if one determines the set of osculating elements at 
one time, say T, and then uses that set as described above to advance 
the position of the spacecraft along the orbit, then considerable error 
in the predicted Cartesian state may develop quite rapidly. The reason 
for this is that there may be considerable deviation in the value of a 
given element in different points of the orbit. (Compare points a and 
b in fig. 3-3.) For example, in a spacecraft orbit of small eccentricity 
at an altitude of 600 km and inclination of 57°, the mean value of the 
semimajor axis is 6981 km. However, the actual value can range from a 
minimum value of about 6976 km to a maximum of 6986 km. Table 3-1 
shows the differences in mean anomaly that one would compute at 
various time intervals using these values of a. 
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Figure 3-3. Sketch illustrating various orders of perturbations. 


Table 3-1. Change in Mean Anomaly for ±5-km 
Change in Semimajor Axis 


a, km 

n, deg/min 

M at t = 100 min, deg 

6976 

3.725059 

372.5059 1 

6981 

3.721058 

372.1058 

6986 

3.717064 

371.7064 


Since 1 represents a down track distance of about 122 km for this 
orbit, table 3-1 represents a potential difference of about 97 km in only 
100 min, which is about 1 revolution for this orbit. Even if the first- 
order perturbations are included as described in chapter 4, rather large 
errors still ensue after only a short time. 

Experience has shown that much greater accuracy is maintained 
over long time periods if the mean elements are used. Perhaps a better 
way of saying this is to say that a more acceptable error results which 
allows the simple algorithm to be used over a longer period of time. 
Relations for the osculating elements in terms of the mean elements 
can be found in the references by Kozai (1959), Brouwer (1959), et al. 
cited earlier. These same relations can be used iteratively to compute 
the mean elements from the osculating set at any time. 
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Chapter 4 

Equations of Motion 

The development of the general equations of motion for a satellite 
of negligible mass orbiting about a massive primary, or central body, 
can be found in practically any textbook on Celestial Mechanics. (See, 
for example, Escobal 1965, Smart 1977, Dubyago 1961, or the classical 
text by Moulton (1914).) In vector form 

$ = W H-') 

where r is the radius vector from the center of the primary (the origin of 
coordinates for this development) to the spacecraft, t is time, and V is 
the gravitational potential of the central body. If we restrict ourselves 
to the 6th-degree zonal expansion of equation (2-9) and recall that T~ 

sin V 7 “ 
r 

then we can write the three components of equation (4-1) as (see, for 
example, Escobal 1965) 
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(4-3) 
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(4-4) 


and according to GSFC, the latest (as of March 1986) values of the 
constants are 

M = 398600.64 km 3 /sec ^ 

R e = 6378.140 km 
J 2 = +1082.28^-6 
J 3 = — 2.5358868£-6 
J 4 = —1.6246180^-6 
J 5 = —0.22698599^-6 
J 6 = +0.54518572^-6 

Equations (4-2) to (4-4) can readily be programmed on a modern 
computer and integrated using any number of accurate, rapid numerical 
methods. A modified 7th-order Runge-Kutta technique, one of several 
available canned routines at the Langley Research Center (LaRC), was 
used for the numerical examples of the present text and was also used 
in some of the operational versions of the integration routines. 

The equations of motion can also be written directly in terms of 
the orbital elements (the Lagrange Planetary Equations.) (See for 
46 


I 


J3 

t 




I?- l’ 


t ; t t 


-t -l • U V 


I L l ’* 


r .r 





Chapter 4 

example, Escobal 1965, Smart 1977, Dubyago 1961, Moulton 1914.) 
The perturbative gravity field, however, becomes very messy when 
expressed in terms of the orbital elements. Consequently, this form 
of the equations of motion, while occasionally used in numerical work, 
is mainly used in theoretical developments, in which usually only the 
J 2 or the J3 term of the zonal expansion is used in an attempt to 
derive purely analytical expressions for the response of the spacecraft 
orbit to these terms. (See, for example, Kozai 1959, Brouwer 1959, 
Garfinkel 1959, King-Hele 1958.) Escobal (1965) expounds somewhat 
on Kozai’s method and gives a lucid physical explanation of the effects 
of the perturbations on the orbital elements. (To this end, see also the 
excellent presentation by Moulton (1914).) 

As mentioned earlier, the perturbing potential is difficult to express 
in terms of the orbital elements and, with modern computers and 
numerical techniques, it is the author’s contention that, if one is merely 
seeking accurate numbers with which to work, then the integration of 
equations (4-2) to (4-4), or in their more complete form equation (2-6) 
for the potential, is as good a way as any. Of course, if one is faced 
with the problem of predicting the position or orbital characteristics 
far in the future — for example, for several weeks or months — obviously 
integration of the equations of motion would become quite expensive. In 
this case, one would probably be well advised to sacrifice some accuracy 
for economy and speed and use one of the theoretical models alluded to 
earlier or the simple mean element model shown later in this chapter 
which uses only the J 2 term. 

The first numerical examples are shown to illustrate the accuracy 
of equations (4-2) to (4-4). An initial state vector 

x = 3211.365 km 
y — -4680.423 km 
2 = —4081.154 km 
x = 2.326315 km/sec 
y = 5.555629 km/sec 
z = —4.545389 km/sec 

was picked from an ephemeris tape prepared by GSFC for the SAGE II 
experiment. The GSFC prediction program, as mentioned earlier, uses 
an n = m = 8 gravity field (eq. (2-6)) and serves as a standard here 
with which other computations are compared. 
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The same initial conditions were used in the integration of equa- 
tions (4-2) to (4-4) with three runs: 

1. All J 2 -J& used 
2- J 2 ± 0, J 3 ~J e = 0 
3. All J 2 ~J 6 = 0 

Run 3, of course, corresponds to the purely central force field. The 
results of these runs are given in table 4-1 as follows: the first column 
lists the GSFC results after 48 hr (2880 min), 96 hr (5760 min), and 
144 hr (8640 min). Columns 2-4 list the differences 

(LaRC results of runs) - (GSFC results) 

and the LaRC results were computed from the modified Runge-Kutta 
numerical method with a 60-sec time increment for the computation 
interval. 


Table 4-1. Comparison Between GSFC and LaRC Ephemerides 

T Run 1 includes all J 2 -J 6 ; Run 2 uses J 2 only; ] 
l Run 3 is spherical Earth result 



GSFC 

Differences in LaRC and GSFC results for— | 

| Run 1 

1 Run 2 

Run 3 


t = 48 hr 


x, km . . . 

-2414.451 

-3.087 

-3.686 

591.649 

y, km . . . 

-5520.263 

4.505 

5.144 

-57.877 

z, km . . . 

3521.274 

4.912 

4.763 

257.880 

x, km/sec 

3.1777850 

-0.002647 

-0.002981 

0.504577 

y, km/sec 

-4.6332890 

-0.006146 

-0.006592 

0.140997 

i, km/sec 

-5.0583560 

0.003965 

-0.003546 

0.223293 


t = 96 hr 


x, km . . . 

-2767.378 

4.915 

5.892 

-788.005 

y, km . . . 

3806.603 

10.928 

12.370 

-212.267 

z , km . . . 

5141.751 

-4.985 

-4.834 

-341.438 

i, km/sec 

-2.997119 

-0.005224 

-0.006453 

1.376789 

y, km /sec 

-6.258721 

0.007867 

0.008922 

-0.183703 

i, km/sec 

3.014856 

0.010240 

0.009824 

0.608086 

t = 144 hr 


x, km . . . 

3164.478 

6.164 

7.876 

-2011.458 

y, km . . . 

5901.433 

-7.421 

-6.162 

354.769 

x, km . . . 

-1980.466 

-13.117 

-11.255 

-894.553 

x, km/sec 

-2.952138 

0.008289 

0.009102 

-1.038673 

y, km/sec 

3.573608 

0.015006 

0.013410 

-0.310911 

z y km/sec 

5.968511 

-0.005320 

-0.005758 

-0.443209 
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It is seen from table 4-1 that the use of only the J2~J% terms of the 
gravity field expansion is not too shabby when compared with the full- 
blown 8 by 8 model of GSFC and is certainly accurate enough for many 
satellite applications. Equally obvious from the table is the fact that 
this model is probably not accurate for other applications — for example, 
long-range prediction (i.e., several months) of the ground or geodetic 
position of the spacecraft to be used for ground truth location studies. 
However, even for 6 days of integration, the differences between the 
standard GSFC model and the truncated model used herein are only 
of the order of 10 km or so in position and of the order of 15 m/sec in 
velocity. 

The J 2 term clearly dominates the perturbation effect, a fact which 
has been known since the early days of satellite flight. In fact, 
examination of the table shows that, at least over the 6-day period 
used for the example, the use of the J$-J§ coefficients is not warranted 
at all. However, only further study could determine the time span in 
which their exclusion should be adopted; for example, their effect would 
doubtless show up for an integration span of several weeks. 

If one now converts the Cartesian components of table 4-1 to 
orbital elements using the CARCON algorithm, one finds the results 
of table 4-2. Here, the complete values are presented, not merely the 
differences between the LaRC and the GSFC results. The actual orbital 
elements at t = 0 hr were 


a = 6981.471516 km 
e = 0.00141817 
i = 57°. 002219 
0 = 96°. 623064 
u> = 58°. 316978 
M = 165°. 753617 

First, note the last column of table 4-2, which shows the unper- 
turbed elements. The agreement of these elements (except, of course, 
M) with the initial conditions illustrates the accuracy and stability of 
the numerical integration routine and serves to rule out integration er- 
rors as a cause of some of the differences in the Cartesian coordinates 
of table 4-1. 

Second, the orbital elements appear to be very close in all cases. 
The percent differences in all elements except u ; and the anomalies are 
very small (order of 0.001 percent) compared with the percent error 
in the Cartesian coordinates (order of 0.2 percent). The maximum 
errors are in and / individually. However, and much more important, 
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the errors in the sum u + / are rather small, and it is the sum 
which determines the angular distance of the spacecraft from the 
ascending node. The error in u ; is mainly due to the small eccentricity 
of this orbit. As the CARCON algorithm of chapter 3 shows, the 
computation of u = cu + / does not explicitly involve the eccentricity, 
but the computation of / does. Therefore, u ought to be determined 
with reasonable accuracy, even for small eccentricities. However, when 
the eccentricity is small in absolute value, even very small absolute 
errors are very large relatively, or percentagewise, and hence the 
accuracy to which / , and consequently u, are individually computed is 
questionable. If these are used to determine the Cartesian coordinates, 
though, the results are still reasonable, as u is the quantity which is 
used, and this is determined with some confidence. 

If one considers only the J 2 part of the perturbed gravity field and 
further restricts oneself to first-order theory, then one finds that a, e, 
and * are constants in time for a given orbit. The longitude of the 
ascending node, ft, and the argument of perigee, w, are found to vary 
linearly and are given by the following equations: 

n — n Q + n (t to) (4-5) 

= UJo + u(t - to) (4.6) 

M = M 0 + n(t — t 0 ) (4-7) 

where 



and in which 

p = a( 1-e 2 ) (4-11) 

If equations (4-5) and (4-6) are used to compute fi and uj and equa- 
tions (4-7) and (4-10) are used to compute M, it is found that the 
element-to-Cartesian coordinate algorithm given earlier will predict the 
Cartesian components reasonably accurately for longer periods of time 
than can be computed from the central force alone (order of a few 
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days), but of course the acceptable time period depends heavily on the 
intended use and accuracy requirements. 

Table 4-2. Comparison Between GSFC and LaRC Osculating 
Orbital Elements Computed at 48, 96, and 144 hr 


[ Run 1 includes all Run 2 uses J 2 only; 

Run 3 is spherical Earth result 



GSFC 

LaRC results for — 

Run 1 

Run 2 

Run 3 

t = 48 hr 

a, km 

6983.088680 

6983.070167 

6983.090200 

6981.471471 

e 

0.00186271 

0.00184027 

0.00181924 

0.00141816 

»\ deg 

57.006807 

57.006451 

57.006509 

57.002221 

0, deg .... 

88.674902 

88.675025 

88.667845 

96.623060 

w, deg .... 

72.635991 

72.802224 

74.944251 

58.318410 

M, deg .... 

70.179204 

69.955204 

67.817778 

81.313131 

f, deg .... 

70.380168 

70.153465 

68.010984 

81.473819 

w + /, deg . . . 

143.016159 

142.955684 

142.955235 

139.992229 

t = 96 hr 

a, km 

6977.549839 

6977.616167 

6977.629441 

6981.472607 

e 

0.00116235 

0.00113651 

0.001140987 

0.00141831 

*, deg 

56.992038 

56.991955 

56.991988 

57.002216 

Q, deg .... 

80.804759 

80.806862 

80.792607 

96.623066 

w, deg .... 

89.510732 

91.289731 

98.225995 

58.319535 

M, deg .... 

332.149893 

330.265486 

323.337990 

356.872952 

/, deg .... 

332.087589 

330.200812 

323.259833 

356.864070 

w + /, deg . . . 

61.598321 

61.490543 

61.485828 

55.183605 

t - 144 hr 

a, km 

6986.459458 

6986.385047 

6986.376400 

6981.471158 

e 

0.00146659 

0.00140939 

0.00127645 

0.00141834 

deg 

57.014001 

57.015303 

57.015294 

57.002222 

f2, deg .... 

72.866127 

72.869378 

72.847804 

96.623062 

w, deg .... 

50.919996 

52.026907 

58.879086 

58.318779 

JVf, deg .... 

289.476609 

288.229329 

281.390348 

272.434666 

/, deg .... 

289.318070 

288.075816 

281.246913 

272.272271 

w + /, deg . . . 

340.238066 

340.102753 

340.125999 

330.591050 


Equations (4-5) to (4-10) were used to generate the next set of 
numerical data, again for t = 48, 96, and 144 hr. The set of orbital 
elements given above was assumed at t = 0. The orbital elements 
were updated as prescribed, and the CONCAR algorithm was used 
to generate the Cartesian components from the orbital elements. The 
results are shown in table 4-3. The first two thirds of the table give the 
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results of the appropriate calculations, the top third giving the orbital 
elements and the middle third giving the Cartesian coordinates. The 
bottom third gives the difference in the Cartesian coordinates between 
the above results and the reference GSFC results from column 1 of 
table 4-1. The bottom third also displays the difference between the 
sum w 4- / computed from the exact Cartesian coordinates and that 
computed from the approximate orbital elements. 

The differences shown in table 4-3 are about an order of magnitude 
greater than those found using only J 2 with the integration routine. 
(The integration of the equation of motion includes all orders of effect 
and not just the first-order effects predicted by equations (4-5) to (4- 
7), with a, e, and i constant.) Comparison between tables 4-2 and 4-3 
shows that the nodal point Q is predicted rather well by equation (4-5), 
and the quantity u + f is the prime error source. However, we 
must remember that the approximate orbital element treatment using 
equations (4-5) to (4-10) is expressly to be used with mean orbital 
elements, while we are here applying these equations to osculating 
elements. There may be some instances where we have no choice as 
the only initial conditions available may be in Cartesian form, and 
these almost always refer to osculating elements. Proceeding with this 
caveat, note that the two orbital planes line up rather well but the 
position of the spacecraft in the true orbit is slightly ahead of the 
position in the approximate orbit. This is, in turn, traced to the fact 
that the first-order theory assumes that a is a constant. However, for 
this orbit, a varies between 6975.2 and 6987.6 km, and hence, the tilling 
between the true and approximate (through the mean anomaly) orbit 
is off somewhat and, as can be seen in the last column of table 4-3, 
the true spacecraft position is gaining on the approximate position bv 
about 0.25 deg/day. 

The approximate first-order analysis can thus also be reasonably 
accurate for a few days or so. If one starts off with a set of mean 
elements, then the spatial position can be predicted reasonably well for 
periods of several weeks to several months, depending on the accuracy 
constraints imposed on the desired results. This algorithm is rather 
simple to program and compute when compared with the integration 
routine described earlier. 

Table 4-4 demonstrates the application of this modified algorithm. 
The mean elements at t = 0 were computed using the method of 
Kozai (1959). (These equations are rather lengthy and hence are not 
reproduced here. The original reference should be consulted.) The 
mean elements a , e, and i were then held constant, and the remaining 
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three mean elements updated using the following rates computed from 
equations (4-5) to (4-11): 

n = 223.234095 deg/hr 
fl = -0.16475043 deg/hr 
u = 0.073098627 deg/hr 

Table 4-3. Orbital Elements, Cartesian Coordinates, and Errors in 
Cartesian Coordinates Computed With CONCAR Algorithm Using 
J 2 Term and Starting Conditions as Osculating Elements at t = 0 



t = 48 hr 

t = 96 hr 

t = 144 hr 

a, km 

6981.471516 

6981.471516 

6981.471516 

€ 

0.00141817 

0.00141817 

0.00141817 

t\ deg 

57.002219 

57.002219 

57.002219 

ft, deg 

88.712177 

80.801189 

72.890102 

w, deg 

61.826844 

65.336711 

68.846577 

M, deg 

80.515297 

355.276977 

170.038657 

w + /, deg 

142.502536 

60.600278 

338.722663 

x y km 

-2437.812650 

-2718.177361 

3232.649280 

y, km 

-5484.269655 

3907.133828 

5811.739703 

z t km 

3563.455210 

5094.024327 

-2124.784121 

i, km/sec 

3.157658 

-3.050433 

-2.863125 

y, km/sec 

—4.681182 

-6.184601 

3.739829 

z, km/sec 

-5.023572 

3.114084 

5.908481 

Ax, km 

23.362 

—49.201 

-68.171 

Ay, km 

-35.993 

-100.201 

89.693 

Ax, km 

-42.181 

-47.717 

144.318 

Ax, km/sec .... 

0.020127 

0.053314 

-0.89013 

Ay, km/sec .... 

0.047892 

-0.074120 

-0.166221 

Ax, km/sec .... 

-0.034784 

-0.099828 

0.600030 

A(u; + /), deg . . . 

-0.513623 

-0.998043 

-1.515403 


V 


\ 


T 


The top third of table 4-4 gives the mean elements at t = 0, 48, 96, and 
144 hr. The middle third of the table gives the Cartesian coordinates 
computed from the CONCAR algorithm and using the mean elements, 
and the bottom third gives the differences between the LaRC and the 
GSFC results as was done in tables 4-1 and 4-3. It can be seen that 
these differences in table 4-4, using the mean elements at t = 0, are not 
substantially different from those of runs 1 and 2 of table 4-1, being 
within a factor of 2 or 3 of the integrated results, and certainly much 
better than the 1 to 2 orders of magnitude differences displayed in 
table 4-3 in which the oscillating elements were used at t = 0. It 
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can perhaps be implied that the perturbation-modified algorithm using 
mean elements can be used, with only a slight degradation in accuracy, 
over approximately the same span as the numerical integration using 
only the zonal harmonics. It may just be a personal prejudice, but 
the author believes that the numerical integration algorithm is more 
defensible than the modified mean-element approach, especially if the 
integration time exceeds a few weeks. However, the mean-deviant 
algorithm has been successfully applied to the prediction of both 
SAGE II and SAM II ground truth sites with lead times of as much 
as 4-6 weeks. 


Table 4-4. Orbital Elements, Cartesian Coordinates, and Errors in 
Cartesian Coordinates Computed With CONCAR Algorithm Using 
J 2 Term and Starting Conditions as Mean Elements at t = 0 



t = 0 hr 

t = 48 hr 

t = 96 hr 

t = 144 hr 

a, km 

6981.26555 

6981.26555 

6981.26555 

6981.26555 

e 

0.00254626 

0.00254256 

0.00254626 

0.00254626 

deg 

56.997801 

56.997801 

56.997801 

56.997801 

H, deg 

96.601960 

88.693939 

80.785919 

72.877898 

deg 

71.220024 

74.728758 

78.237492 

81.746226 

M, deg 

152.821231 

68.057791 

343.294351 

258.530911 

x, km 

2215.11242 

-2409.67576 

-2755.911096 

3176.32012 

km 

-4679.87474 

-5520.91264 

3819.199893 

5889.11718 

2 , km 

-4089.14043 

3515.54137 

5130.239109 

-2004.32760 

x, km/sec 

2.32563390 

3.18182136 

-3.00963556 

-2.94169646 

y, km/sec 

5.55145098 

—4.63367603 

-6.25485145 

2.59164655 

i, km/sec 

-4.53986239 

-5.06051117 

3.03212667 

5.95682454 

Ax, km 

3.747 

4.775 

11.467 

11.842 

Ay, km 

0.548 

0.133 

12.597 

-12.316 

A 2 , km 

-7.986 

-5.733 

-11.512 

-23.868 

Ax, km/sec .... 

-0.0006811 

0.0040364 

-0.0125166 

0.0104415 

Ay, km/sec .... 

-0.004178 

-0.003870 

0.0038695 

0.0180386 

A 2 , km/sec .... 

0.0055266 

-0.0021552 

0.0172707 

-0.0116862 
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Chapter 5 
Where Is the Sun? 

Most of the time series relations needed to determine the position 
of the Sun have already been given in chapter 1. The remaining 
geometrical concepts will be developed in this chapter, as needed. 

Perhaps the best way to introduce this material is through numerical 
examples. The following problems, typical of those encountered in 
spherical astronomical applications, will be addressed in this chapter, 
and a complete numerical example will be worked out for each to 
illustrate the method. The text procedures are not the only approach, 
of course, but are the ones that, for now unknown reasons, appealed 
to the writer at the time he first encountered these problems. The 
number of significant figures shown in the numerical results is generally 
that required for the accuracy given; 0.1 arc-sec is 0° .000028, and 
calculations to a few hundredths of an arc-sec are frequently required. 

The most fundamental problem which arises concerning any compu- 
tation of the Sun’s position is the determination of the right ascension 
and declination of the center of the Sun, given any date of the year and 
time of day. These coordinates, along with the Greenwich sidereal time, 
are central to the solution of all the problems listed below. Therefore, 
the problems to be considered in the present chapter are the following 
(all the problems presented require the day of the year and/ or the time 
of day to be given; some also require the specification of a particular 
location of the Earth’s surface — the coordinates of an “observer;” these 
are assumed given): 

1. Given any calendar date in the year and time of day, what 

are the right ascension and declination of the Sun? Note: Only 
mean values will be determined as the extra computation will 
serve no useful purpose here. 

2. What is the Greenwich sidereal time for a given universal time, 
and what is the local sidereal time at a specified geographic 
location? 

3. What are the geographic coordinates of the subsolar point at a 
given time of day? 
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4. For a given observer, what are the azimuth and elevation 
angles of the Sun at a given local (zone) time of day (with and 
without refraction)? 

5. What is the local zone time when the Sun is directly on the 
local meridian (zonal time of local noon)? 

6. What is the elevation angle of the Sun at local noon? 

7. What are the local zone times of sunrise and sunset? 

8. What are the azimuth angles of the Sun at sunrise and 
sunset? 

Note that if the Cartesian coordinates of a spacecraft are known in 
the inertial coordinate system, then its right ascension and declination 
can be found from equation (1-la). Then problems 2-8 can be applied 
to satellite motion by substituting the satellite coordinates for those 
of the Sun, taking into account, when necessary, the much more rapid 
time changes of these coordinates. 

As a specific application of some of the concepts discussed here — 
given the semimajor axis and eccentricity of a satellite orbit, what must 
be the location of the line of nodes of a Sun-synchronous orbit which 
crosses the equator at a specific local time? (See, e.g., Brooks 1977.) 

The numerical example will be worked out for the time of April 6, 
1985, at 2:37 PM EST in Hampton, Virginia, whose latitude and 
longitude, respectively, are approximately 37°N and 76°W. 

The chapter will terminate with the solution of a problem arising 
in connection with Sun-viewing satellites (e.g., SAM, SAGE) and 
specifically the problem of reducing the data from these satellites— 
namely, what is the minimum altitude of a ray from the center of the 
Sun (or indeed any other specific celestial location) to the spacecraft 
above the oblate Earth, and what are the geographic coordinates of the 
“subtangent” point on the surface of the Earth? (See Brooks 1980.) 

Problems and Examples 

Problem 1: 

Figure 5-1 (a) shows the inertial coordinates defined earlier in chap- 
ter 1. The coordinates of the Sun are the right ascension a and decli- 
nation 6. Figure 5-l(b) shows the spherical triangle -ySQ. 

For the given date, find the Julian date as outlined in chapter 1. 
Then proceed by the following steps: 

(a) Compute the mean longitude of the Sun, L = M + Cb, from 
equation (1-19) 
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(a) Sun in inertial coordinate system and its relation to Greenwich meridian. 



(b) Blowup of spherical triangle 6SQ of figure 5- 1(a). 

Figure 5-1. Location of Sun in both inertial and Earth- fixed coordinate systems. 

(b) Find the mean anomaly of the Sun, Af, from equation (1-39) 
and the eccentricity of the solar orbit from equation (1-38) 

(c) Compute the quantity J - M from equation (3-7) and convert 
to degrees 

(d) Then, (w + /) = (w + Af) + (/ - Af) 

(e) Compute the obliquity of the ecliptic, £, from equation (1-21) 
Then from figure 5- 1(b), the mean right ascension is 

tan a = cos e tan (Cj + /) (5-1) 

and the mean declination 

sin <5 = sin e sin (Q + /) (5-2) 
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The quadrant of a is always the quadrant of (u> + /). 

Example: 

(l-ll) 3 JDO = 728495 - 3475 + 122 + 6 4- 1721013.5 

= 2446161.5 at 0 hr GMT on April 6, 1985 

Hampton, Virginia, is 5 time zones from Greenwich — <76/15> — and 
therefore 2:37 PM EST is 2:37 + 12 + 5 = 19:37 GMT. Then, 

JD = 2446161.5 + 19.6167/24 = 2446162.31761 
(M3) T u = 0.852630181 

Note on accuracy. T u can be accurately computed from equation (M3) 
if one groups the terms as 

T = ((JDO - 2415020.) + UT/ 24] 

“ 36525. 

Continuing 

(1-19) L = 15°.0390181 

(1-39) M = 92°. 35203707 

(1-38) e = 0.0167156694 

(1-21) £ = 23° .44119896 

(3-7) f -M = 1°.911865208 

(Step (d)) / + u = 16°. 95088331 

(5-1) a = 15°. 62304219 = l*02 m 29 a .5 

(5-2) 6 = 6°. 660242901 = 6°39'36".87 

The procedure described here gives the right ascension and dec- 
lination relative to the mean equator and equinox of date. The ap- 
parent right ascension and declination (referred to the true equator 
and equinox of date) may be obtained from the mean values by ap- 
plying the proper corrections for nutation and planetary aberration. 
(See, for example, Smart 1977, chaps. 8 and 10.) The numerical values 
tabulated in the Astronomical Almanac are apparent values. Linear 
interpolation in AA85 gives a = l /, 02 m 27 s .5 and 6 = 6°39'25".56; 
thus, for this example, Aq = (29.5 - 27.5) x 15 = 30 arc-sec, and 
AS — 36.87 — 25.36 = 11.51 arc-sec. 

Planetary aberration (correction for the finite speed of light) gives 
a correction of about —20" ± 2" on the right ascension, where the 

The number in parentheses at the left of the equation is the number of the 
equation used to calculate the quantity. 
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correction of 2 /r is sinusoidal with a period of 1 year. The aberration 
correction to declination is also a yearly sinusoid with an amplitude of 
about 8 n . The nutation in right ascension is a sinusoid with a period 
of just over 18 yr with an amplitude of about 20 ,f . The nutation in 
declination is a multiple-frequency sinusoidal sum, but its maximum 
magnitude is about 9". 5. Using the formal equations (e.g., Smart 1977, 
chap. 8), these data for the present example are 

Nutation Aberration Sum 

Aa . . . —12". 74 —18". 92 -3l".66 

AS . . . —4". 06 —8". 02 -12". 08 

If these are added to the mean values computed above, we get 

a = l**02 m 27®39 
6 = 6° 39' 24". 29 

which are very close to the tabulated AA85 values. 

For most Earth orbital applications, it is probably not necessary to 
add the nutation and aberration corrections; the mean values should 
suffice. j ~ 

Problem 2 : 

Recall that the sidereal time is the hour angle of the vernal equinox. 

If it is measured from the Greenwich meridian, it is known as Greenwich 
sidereal time, and if it is measured from any other location, it is known 
as local sidereal time. Figure 5*2 shows the Greenwich meridian, with 
the Greenwich sidereal time 6 g . The observer at O has an east longitude 
A 0 and latitude ip 0 . 

(a) Find the sidereal time at 0 hr GMT from equation (1-12) 

or (1-15) using the appropriate T u 

(b) Correct for time of day using equation (l-14a,b) 

(c) From figure 5-2, the local sidereal time is 

LST = 6 g + A c (A 0 + East) (5-3) 


Example: 

JDO = 2446161.5 
(M3) T u = 0.8526078029 

(1-12) 0g o = 194°. 2277554 = 12 h 56 m 54 a .6613 (AA85 gets 
12 fc 56 m 54*.7273) 

= 19 h 37 m GMT = 1177 min 
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(l-14a) A 6 g = 295°. 0556212 

(l-14b) 6 g = 129°. 2833766 

(5-3) LST = 53°. 28337659 = 3 /l 33 m 08 s .0100 



Figure 5-2. Diagram similar to figure 5-1 but including observing meridian. 
Defines azimuth A and zenith distances z of celestial body S relations to 
observer O. 

Problem 3: 

From figure 5-2, the east longitude of the subsolar point, A a , is found 
from 

a = 6 g + \ s ( 5 . 4 ) 

The geocentric latitude of the subsolar point is identical to the decli- 
nation, and the geographic latitude is found from equation (1-3) with 
o - 6378.160 km and b — 6356.775 km. 

Example: 

(5-4) X s = 246°. 3396656 

= 6° .660242901 
(1-3) tpg = 6°. 704796079 

As a rough check, note that the Sun is (360° - 76°) - 246° .3 or about 
38° west of Hampton at this time. If we assume that the Sun is directly 
overhead at 1200 hr local time, which is not exactly true, then in 2 / *37 m 
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the Sun will move to the west, at 15 deg/hr, a distance of about 39°. 25, 
which is close to the 38° above. 

Problem 4: 

In figure 5-2 at the observer’s meridian O, the azimuth angle A is 
measured clockwise from geographic north. The solar zenith angle z , 
arc OS in the figure, is the complement of the altitude angle, and hence, 

7 = 90 ° - * (5-5) 

The zenith angle z is found by applying the law of cosines for sides to 
the spherical triangle POS (the quadrant of z is no problem, since if 
cos z < 0., the object (Sun) is below the horizon): 

cos z — sin i p 0 sin 6 -I- cos cos <5 cos (A s — A 0 ) (5-6) 

Both sine and cosine of the azimuth angle must be used to determine 
the quadrant of A from the triangle, 

cos 6 sin (A s — A 0 ) 

sin A = : 

sin z 

and also 

sin 6 — cos z sin 

COS A = : 

sin z cos tp 0 

Note: Some users prefer to measure azimuth east or west of due south — 
most spherical astronomy texts define it this way. This azimuth A f is 
found from the above as 


(5-7) 

(5-8) 


A' = A - 180° (5-9) 

and is west if A ' is positive, east if A' is negative. 

Example: 

(5-6) z = 45°. 7516467 

(5-5) A = 44°. 2483533 

(5-7) sin A = -0.8471830565 

(5-8) cos A = -0.5313011094 

A = 237°. 9065922 
(5-9) A' = 57° .90659223 W 

To correct for the effects of refraction, the following law, derived for 
a spherical atmosphere (Smart 1977, Green 1985, for example), 


o 

2 — Zft = R = A tan zr + B tan z^ 


(5-10) 
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is used, where the constants A and B, while defined by analytic 
expressions, are usually found empirically. Smart (1977) gives 

A = 58". 16 

B = — 0".067 


whereas Green (1985) gives 

A = 60". 29 


B = -0". 06688 


The above constants assume standard atmospheric conditions at the 
surface. Green (1985) gives the following expressions for A and B, 
in which the temperature and pressure profiles of the atmosphere are 
allowed to vary. Define 

roc 

PoHo = I p dh 

Jo 


where p is the density at altitude h. The parameter H 0 is thus the 
density scale height of the atmosphere. Then Green gives the correction 
equations 


A= (n Q - 1) 



(5-11) 


B = — (n 0 - 1) 


Ho 1 , 


(5-12) 


where the units of A and B are radians, and n 0 is the index of refraction 
of sir at the surface. This can be computed with acceptable accuracy 
from the well-known Edlen formula (Edlen 1953) 

r>6 _ („ , 0.459 \ P 

A 2 ) T 

(„ « _ (5 ' 13> 


(no - 1) x 10 6 = ^77.46 + 


1013 


43.49 - 


A 2 ) 


where 

P = atmospheric pressure, mb 
T = atmospheric temperature, K 
Ph 2 0 = P art ial pressure of water vapor, mb 
A = wavelength, pm 


Given R, then the observed zenith as seen at the surface of the Earth 
is 
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It is seen that if extreme accuracy is required, an iteration is required. 
Assume zr = then compute z from equation (5-6). Use this 
in equation (5-10) to get an estimate of R\ if further accuracy is 
needed (seldom), use zr computed from equation (5-14) to go back 
to equation (5-10) and recompute R. 

The refraction correction described above (equally applicable, of 
course, to the moon, stars, satellites, or any other celestial body) is 
sufficiently accurate whenever zr < 75°. For zenith angles greater than 
this, tables of refraction correction are generally used. The following 
equation (see any Astronomical Almanac, Section B) has been quoted 
as being reliable for zr > 75°. Let 7 equal 90 — z , the altitude. Then, 

(° 1594 + 0.01967 + 0.000027 2 ) 

” (273 + T)(1 + 0.5057 + 0.08457 2 ) ’ 

where 


P = surface pressure, mb 
T = surface temperature, K 
7 = altitude, deg 

and in this formula, R is in degrees. For z = 75°, and at STP 
(P = 1013 mb, T = 0°), the above formula gives R = 0°. 06159 or 
221". 724. Equation (5-10) gives 213". 573 using Smart’s constants and 
221". 529 using Green’s. 

Example: Use Smart’s constants. 

Assume zr = 45.7516467 
then 

(5-10) R = 5 9". 633854 
(5-14) z R = 45°. 735081 = 45°44'06".29 
The new R would be 59". 599454 and the iterated zr = 45°. 735091 or 
45°44 / 06".33, or a change in zenith distance of 0".04. 

Problem 5: 

In figure 5-2, the observer’s meridian is rotating from west to east 
at the rate of 15 deg/mean sidereal hr. From the figure it is obvious 
that the condition 

Qg + Ao — QL 
or 

®9o + ^min + A 0 = Q (5-16) 

is the condition for the Sun to be on the local meridian. We know that, 
within a maximum error of 7°. 5 or 30 min in time, the mean Sun is 
directly overhead at 12 hr LMT. Thus for EST, for example, local noon 
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will occur near 17 hr GMT. Let us first find the right ascension of the 
Sun at 17 hr GMT and assume that it is fixed at this position for this 
calculation. Proceeding exactly as in the first example, for 17 hr GMT, 

a = 15°. 52331572 
6 = 6°.619146553 

From equation (5-16), 

194°. 2277554 + 0.25068447f min + 284°.0 = 15°.52331572 

Solve for f min (using arithmetic mod 360), we find that t min = 
1026.3725 min after 0 hr GMT, or 17 /l 06 m 22 s GMT. Since our ob- 
server is in an Eastern Standard Time zone, the local time of the Sun’s 
meridian passage is 12 A 06 m 22*. Interpolation in the tables given in 
AA85 gives 12 / *06 m 29 s , and our approximate method here is only 6 sec 
in error (or else AA85 is 6 sec in error). 

Problem 6: 

From equation (5-6), at local noon, X a = A 0 , and 

cos z n oon = cos %> 0 cos 6 + sin ipo sin 6 = cos {ipo ~ b) 

^noon = *Po — b 


Example: 

(5-11) Znoon = 37° .0 - 6° .660242901 = 30°.3397571 

(5-5) T'noon = 59°. 6602429 
Problem 7: 

The angle measured clockwise (looking down at the north pole) 
between the observer’s meridian and the meridian of any celestial body 
is called the local hour angle of that body — for example, in figure 5-2, 
the angle OP7 is the local hour angle of the vernal equinox, H~, and 
the angle OPS is (360 — H s ), where H s is the hour angle of the Sun; 
that is, 


360 - H s = A a - A 0 (5-17) 

in figure 5-2. In terms of the hour angle of the Sun, we can write 
equation (5-6) as 

cos z = sin ipo sin 6 + cos rp 0 cos 6 cos H s (5-18) 

If the Earth had no atmosphere, the z = 0° would put the center of 
the Sun right on the horizon, and thus from figure 5-2 
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cos H r j s = - tan \p Q tan <5 (5-19) 

would give the rise/set value of the hour angle, from which the rise/set 
times from the time of local noon could be determined. 

The Earth does have an atmosphere, however, and there is refraction 
to be considered. Consequently, there are three different ranges of 
twilight generally computed in the almanacs: 

(a) Civil twilight begins when the top of the Sun just touches the 
horizon; this requires z = 90°50 / (i.e., 16' for the half-angle of the Sun 
plus 34' refraction correction) 

(b) Nautical twilight begins when z = 96° 

(c) Astronomical twilight begins when z = 102° and ends when 
z = 108° 

To calculate the various times, put the appropriate values of z into 
equation (5-12) to compute the hour angle, divide the hour angle by 15 
to get the number of hours from local noon, and add or subtract this 
time interval from the local time of meridian passage (problem 5) to 
get the respective time of sunrise/sunset. 

Example: 

We know that sunrise/sunset will occur at about 6 hr before and 
6 hr after local noon, respectively. For the EST time zone, these will 
occur at about 11 hr GMT and 23 hr GMT. These times are used as in 
problem 1 to get the declinations of the sun at sunrise/sunset, 

<5 r ise = 6°. 524818313 
fi 8et - 6°. 713346807 

These are used in equation (5-18) to give the values in table 5-1. 
(The numbers in parentheses are the minutes determined by a double 
interpolation in AA85 for the example latitude and times.) 

Problem 8 : 

From equations (5-7) and (5-8), we can now set z = 90°, since here 
we’re concerned with the Sun center and refraction can be neglected, 


sin A = cos 6 sin H 

sin <5 - sin xp 0 

cos A = : 

cos xpo 


(5-20) 

(5-21) 


in which, of course, the appropriate rise/set values are used for H and 6. 


T 
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Table 5-1. Hour Angles and Rise/Set Times Measured 
From the Time of Local Meridian Passage for Example 7 


z 

AM 

PM 

^riae> deg 

^rise> hr 

^set> deg 

*aet> hr 

90° 50' 

95.99935 

6.38248 

96.144453 

6.39213 

96° 

102.58716 

6.82047 

102.73768 

6.83048 

102° 

110.37838 

7.33847 

110.53816 

7.34909 

O 

00 

o 

118.40111 

7.87186 

118.57461 

7.88339 


z 

Sunrise 

Sunset 

90°50' 

5 #i 43 m * (38 m ) 

18^30 m *(28 m ) 

96° 

5 h 17 m (ll m ) 

18 , *56 m (53 m ) 

102° 

4 h 46 m (41 m ) 

19^21 m (24 m ) 

108° 

4 h 14 m (10 m ) 

19 h 59 m (5 7 m ) 


’"Numbers in parentheses are the minutes determined by double interpolation in 
AA85 for the example latitude and times. The present technique appears to be accurate 
to within 5 or 6 min when compared with AA85. 

Example: 

Sunrise 

6 = 6°. 524818313 
H = 94°. 9443029 
sin A = 0.9898252808 
cos A = -0.6112693926 
A = 127°. 681345 

Sunset 

S = 6°. 713346807 
H = -95° .08883614 
A = 232°.6143557 

As one example for the application of this material to an orbit 
problem, consider the following frequently posed problem: a satellite 
is to be placed in a circular (e = 0) orbit at some given altitude (a is 
thus known), such that it is Sun synchronous and the ascending node 
crosses the equator at some specified local time, t hr, either before or 
after local noon. What must be the inclination of the orbit and where 
should the ascending node be placed at orbit insertion? 

Recall that the mean daily angular travel of the Sun about the 
Earth is 360/365.2422, or 0.98564733 deg/mean solar day, and that this 
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motion is west to east in the inertial coordinate system. Equation (4-8) 
gives the angular rate of the ascending node of the orbit: H must thus 
be +0. 98564733 deg/mean solar day in order for the orbit to be Sun 
synchronous. Equation (4-8) then defines the inclination of the orbit, 
which must be greater than 90° in this case — i.e., retrograde — in order 
to be Sun synchronous. (See Brooks 1977, for example.) 

Finally, for a PM crossing, the node must be located to the east 
of the solar meridian (15 x where t is the time from local noon in 
hours), whereas for an AM crossing, the node must be located the same 
angular distance to the west of the solar meridian. 

For instance, suppose we want to put a satellite into a circular Sun- 
synchronous orbit with a semimajor axis of 6978 km. We want to have 
the ascending node cross the equator at 2 PM local time. What are 
the inclination and the right ascension of the ascending node at orbit 
insertion for this case? Assume the insertion date is the same date as 
we’ve used so far for numerical examples. 

At a first guess, assume that ft — n = yj /i/a 3 in equation (4-8). 

This gives n = 5361.7792 deg/day and, with tl = 0.98564 deg/day, 

equation (4-8) yields i = 97°. 789976. To improve on this slightly, use 

this value of i in equation (4-10) and compute n = 5358.3436 deg/day T 

(this is only a 0.6-percent decrease in mean angular rate). Using this 

value for n in equation (4-8) gives i = 97.795001 deg/day, and obviously 

further iteration is unwarranted, as the inclination at orbit insertion 

cannot be achieved with this accuracy. 

On April 6, 1985, at 0 hr GMT, the right ascension of the Sun is 
15°2161669. For a 2 PM local time of crossing, the ascending node must 
be 2 x 15 or 30° to the east of the Sun’s meridian, or the right ascension 
of the ascending node must be at 45.216°. Due to the nonuniform 
motion of the Sun throughout the year, the 2 PM crossing cannot 
be maintained exactly, but neglecting the effects of the higher order 
perturbations, this orbit is the best that can be done. Just as the 
differences in right ascension between the mean and true Sun are given 
by the equation of time, the same equation will predict the differences 
between the actual nodal crossing time and 2 PM. These differences 
will vary from —14 to +16 min throughout the year. Figure 5-3 shows 
a rough sketch of the resulting orbital geometry. 

Omnipresent /? Angle I 

The {3 angle is a parameter that comes up again and again in mission 
design literature, as many mission parameters are directly related to 
this angle or functionally dependent on it (see, e.g., Buglia 1986). The 
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Figure 5-3. Sketch of geometry of retrograde Sun-synchronous orbit, with 
equatorial crossing time t hr PM. 

(3 angle is defined (see fig. 5-4) as the angle between the Sun vector 
and the orbital plane, and as cam easily be seen, many problems related 
to temperature or cooling, Sun viewing or the prevention of solar energy 
entering the optics of am instrument, and many other such mission 
parameters are closely related to this quantity. 



Figure 5-4. Definition of 0 angle. 


If the spacecraft Cartesian coordinates are given in the inertial 
coordinate system, the unit vector normal to the plane of the orbit 
(along the angular momentum vector) is given by 


N = 


r x f 
r x f| 


(5-22) 
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If e s is a unit vector to the center of the Sun, then 


and thus 


e s 


cos 6 cos a 
cos S sin a 
sin 6 


(5-23) 


sin (3 = JV • e $ (5-24) 

In terms of the elements ft and i of the orbit, the normal to the orbit can 
be written in the alternate form from equation (5-14) (Escobal 1965) 
as 


N = 


sin ft sin i 
— cos ft sin i 
cos i 


(5-25) 


Equations (5-23) and (5-25) then permit equation (5-24) to be written 
in the form 


sin (3 = cos i sin 6 + sin i cos 6 sin (ft — a) (5-26) j- 

As stated above, many mission parameters are directly related to the 
(3 angle. For example, solar panels require that they view the Sun, 
whereas some measurement instruments may require that they never 
view the Sun because direct sunlight may either permanently damage 
the instrument or saturate it so that it may become inoperable for a 
short time. If the orbit insertion date is specified, the position angles 
of the Sun are known. Thus, if the (3 angle is specified by mission 
constraints, the right ascension of the orbit can be determined from 
equation (5-26). 

How Long Is the Solar Day? 

The solar day and its variable length were discussed briefly in 
chapter 1, and figure 1-11 shows a plot of the difference between the 
length of the solar day and 86400 sec, the length of the mean solar day. 

We are now able to do the calculations which led to that figure. 

At local noon of some day, the Sun is at position (T) of figure 5-5 ^ 

on the observer’s meridian. The next day at noon, the Sun is again 
in the observer’s meridian at position (2). Thus, during the solar day, 
the Earth has turned through the angle (c *2 — aq) + 360°, where »2 
and a\ are the respective values of the right ascension of the Sun. The 
angular rotation rate of the Earth is 0.25068447 deg/mean solar min 
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(see chap. 1), or 4.178074622 x 10 3 deg/mean solar sec. Thus, the 
length of the solar day is 


(p- 2 ~ <*i) +360° 
4.178074622 x 10" 3 


mean solar sec 


Z 



Figure 5-5. Location of Sun on observer’s meridian on two consecutive days at 
local noon. 

Example: 

FVom figure 1-11 the longest solar day occurs on about December 23 
and the shortest on about September 17. From the methods discussed 
earlier in this chapter we find that from noon on December 23 to noon on 
December 24, the right ascension of the Sun changes by 1°. 1098751, and 
hence the length of the solar day at this time iB 86429.733 mean solar 
sec. Similarly, from September 17 to 18, the right ascension changes 
by only 0°.8967917, and the length of this solar day is 86378.733 mean 
solar sec. The solar day on December 23 is thus about 51 sec longer 
than the solar day on September 17. 

Minimum Height of Ray Above Oblate Spheroid 

The Earth can, for most space applications, be adequately repre- 
sented by an oblate spheroid, with equatorial radius a and polar radius 
b. Suppose we are given a spacecraft position vector r p , and an arbi- 
trary direction from the spacecraft represented by the unit vector e. 
This vector could be, for example, a unit vector to the center of the 
Sun, to another satellite, or to any other celestial body. The question 
we address here is, then, what is the minimum distance between the 
ray defined by e and the surface of the spheroid, Vin. and what are 
the geographic coordinates of the subtangent point A. (See fig. 5-6.) 
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For a spherical Earth, this problem is almost trivial. (See fig. 5-7.) 
The angle p is given by 


p = cos 


-1 T P 


Tn • e 


(5-27) 


and the distance c from the spacecraft to the normal to e is given by 

c = r p cos (180 — p) 

Hence, the tangent point vector rx is 


and thus /i m j n is simply 


r? = r p + ce 


^min r T 


(5-28) 

(5-29) 


( 5 - 30 ) 


where i? c is the mean radius of the Earth. 

In terms of the right ascension and declination of point A, we can 
write 


xj = rx cos 6 t cos ax 
y T = rx cos 6 x sin ax 
zx = rx sin 6 x 

from which the right ascension ax and declination 6x of the subtangent 
point can readily be found. The longitude of A is then found from (see 
problem 2) 



X A = a T - Qg ( 5 - 32 ) 

where 9 g is the Greenwich sidereal time. 

For an oblate spheroid, the problem is still directly solvable, but 
somewhat more cumbersome in form. The author’s solution is presented 
below. 

The vectors r p and e define a plane. In this plane, define new axes x 
and y, where x is parallel to e and y is normal to e. This plane intersects 
the ellipsoidal surface in a curve given by 


f{x,y) = 0 

and the position of the subtangent point (and hence the minimum 
altitude) is defined by the condition 



71 


v 

iL 


vf vi. \ 


-V • U K T t ' l t ’* 

l. 


* • \ ; 


i '■ a: 


i • 
* 



I A An Ann A 


Compilation of Methods in Orbital Mechanics and Solar Geometry 

as can be seen from figure 5-6. The coordinates of A in the “barred” 
coordinate axis system, x A and y A , are then rotated back into the 
“unbarred” coordinate axes (the original axis system), and hence A is 
located in this system. 

Compute r T in exactly the same way as required by equation (5-29), 
using p from equation (5-27) and c from equation (5-28). Then, 
construct the unit vectors 


e T — 


V 

L*J 


is given and 


(5-33) 


t- - r JL - 

m l 


m2 

V r T 

m 3 _ 

to the x-y plane 



n i 

X €y — 


n 2 



. n 3 


(5-34) 


is 


(5-35) 


Now, any vector r(x, y, z) in the unbarred axis system has components 
in the barred axis system given by 


x = r • e± 
V = rey 
z = r • tz 


( 5 - 36 ) 


or, in matrix form, 


X 


h £2 £3 


X 

y 

— 

m, m 2 m3 


y 

z 


n l n 2 n 3 


z 

L m 


(5-37) 


In the present application, we’re interested in the inverse transformation 
to equation (5-37) 


X 


£\ m 1 n. 


X 

y 

= 

£2 m2 n 2 


y 

z 


£ 3 m 3 n 3 


z 


(5-38) 
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Figure 5-6. Geometry of unrefracted ray from celestial body to spacecraft. 

Reciprocal direction is defined by unit vector e; minimum altitude of ray 
(subtangent point) above oblate spheroid is defined along with its geodetic 

and geocentric latitude, <j> g and <f > c , respectively. — . r— i— * — • sr~- 


y 



Figure 5-7. Geometry of unrefracted ray from celestial body to spacecraft. 
Reciprocal direction is defined by unit vector e; minimum altitude of ray 
(subtangent point) along the extension of geocentric radius OA is defined. 


73 


v, 

jL 


is? 

L 




w 


. \ 7 - \ ’ 


A 


U 1 


ir b 


IL l 




i a n a a n 


Compilation of Methods in Orbital Mechanics and Solar Geometry 
and, since we’re interested only in the x-y plane, z = 0, and 


x = i\x + m\y 
y ~ t 2 X + m 2 j/ 
z = £ 3 x + m$y 


(5-39) 


In the unbarred axis system, the equation of the ellipsoid represent- 
ing the Earth is 


or 


x 2 7/ 2 * 2 

~2 H o To = 1 


* 2 + J 2 +p 2 2 = ,> 2 


(5-40) 


If we substitute equations (5-39) into equation (5-40) we get the 
equation of the curve defined by the intersection of the r p -e plane with 
the ellipsoid; that is, f(x, y) = 0, or 


{i\x -|- miy) 2 + (i 2 x + m 2 y) 2 + ^ {t 3 x + m 3 y ) 2 = a 2 (5-41) 
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The coordinates x A , y A of the subtangent point are defined by the 
condition ^ 0. Differentiate equation (5-41) with respect to i, set 

= 0, and collect terms to get 


I A ( £ 1 + *2 + + yA (*irni + t 2 m 2 + ^-£ 3 m 3 ) = 0 (5-42) 


or 


X A = $2/A 


(5-43) 


where 


$ 


VA, 


__ ~ (^ 1™1 + £ 2 m 2 + 

(<?+4 + $4) (5 " 14) 

Use equation (5-43) to eliminate x ^ in equation (5-42) and solve for 


yA = 
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[(fi$ + m,) 2 + (t 2 $ + m 2 ) 2 -I- ^ (^ 3 $ + m3 )2j 1/2 


(5-45) 
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and then x ^ follows from equation (5-43), and the coordinates of the 
subtangent point A in the barred coordinates are known. 

Reference to figure 5-6 shows that the minimum altitude /i m j n is 
given simply by 


Vin = krl - va ( 5 - 46 ) 

Finally, the coordinates of A in the unbarred coordinate axis system are 
found by putting x\ and y\ (with z — 0) into equation (5-38), giving 


*a 


tl mi 


X A 

J/A 

= 

(-2 m2 


VA 

. 2 A . 


_ ^3 m 3 _ 


m 


The right ascension and declination are then found from equation (5-31) 
and the longitude from equation (5-32). 

Of course, the declination is identical with geocentric latitude, and 
finally the geodetic latitude is found from equation (1-3). 

Final Remarks 

It goes without saying that there are myriad other problems, not 
mentioned herein, that are solvable by the methods and equations of 
the present text. With regard to Earth satellites, for example, these 
include such diverse problems as 

1. Rise/set times of a satellite with respect to a given ground 

station 

2. Range, range-rate, and angular data of a satellite with respect 

to a given ground station 

3. Rise/set times of one satellite with respect to another satellite 

4. Range, range-rate, and angular data of one satellite with 

respect to another 

5. Computation of ground tracks of satellites 

6. Area coverage of satellite borne sensors 

7. Entry/exit with regard to sunlight conditions (Escobal (1965 

and 1968) and Green (1985) discuss these and other problems 

in “Keplermanship;” see also Buglia 1986 and Brooks 1977) 

The planets move in essentially elliptic orbits about the Sun. A 
consistent set of the time- varying orbital elements for all the planets is 
given by Escobal (1968). These orbital elements can be used with the 
present equations to determine the heliocentric position of any planet. 
If the orbit of the Earth is also computed, then a few elementary 
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Compilation of Methods in Orbital Mechanics and Solar Geometry 

coordinate transformations are all that is necessary to yield the right 
ascension and declination of the planet with respect to the standard 
Earth- centered coordinates described in chapters 1 and 4, and hence the 
positions of the planets with respect to either a fixed ground station or 
a satellite orbiting the Earth can be readily computed by the methods 
of chapter 5, with the planet’s position being substituted for that of the 
Sun. 

The Moon is the Earth’s nearest natural celestial body, and except 
for the musings of poets, lovers, and other eccentrics, the motion 
of our nearest neighbor is quite complex due to the rather large 
perturbations imposed on its motion by the oblate Earth and the Sun. 

However, Escobal (1968) presents an algorithm, based on earlier work ~- 

by G.W. Hill and E.W. Brown, which is reported to yield an accuracy 
of about 30 arc- sec. The Moon subtends an angle of about V 2 0 as seen 
from the Earth, or a radial displacement of 900 arc-sec. The Escobal 
algorithm is thus accurate to about V 30 of the angular radius of the 
Moon. This is not too shabby and is certainly close enough to find the 
Moon. 

NASA Langley Research Center 
Hampton, Virginia 23665-5225 

June 29, 1988 J~ 
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